{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 数据预处理"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load and clean-up data\n",
    "from numpy import nan\n",
    "from numpy import isnan\n",
    "from pandas import read_csv\n",
    "from pandas import to_numeric\n",
    " \n",
    "# fill missing values with a value at the same time one day ago\n",
    "def fill_missing(values):\n",
    "\tone_day = 60 * 24\n",
    "\tfor row in range(values.shape[0]):\n",
    "\t\tfor col in range(values.shape[1]):\n",
    "\t\t\tif isnan(values[row, col]):\n",
    "\t\t\t\tvalues[row, col] = values[row - one_day, col]\n",
    " \n",
    "# load all data\n",
    "dataset = read_csv('../input/household_power_consumption.txt', sep=';', header=0, low_memory=False, infer_datetime_format=True, parse_dates={'datetime':[0,1]}, index_col=['datetime'])\n",
    "# mark all missing values\n",
    "dataset.replace('?', nan, inplace=True)\n",
    "# make dataset numeric\n",
    "dataset = dataset.astype('float32')\n",
    "# fill missing\n",
    "fill_missing(dataset.values)\n",
    "# add a column for for the remainder of sub metering\n",
    "values = dataset.values\n",
    "dataset['sub_metering_4'] = (values[:,0] * 1000 / 60) - (values[:,4] + values[:,5] + values[:,6])\n",
    "# save updated dataset\n",
    "dataset.to_csv('household_power_consumption.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 模型评估\n",
    "在本节中，我们将考虑如何开发和评估家庭用电数据集的预测模型。\n",
    "\n",
    "本节分为四个部分。他们是：\n",
    "\n",
    "- 问题框架\n",
    "- 评估指标\n",
    "- 训练和测试集\n",
    "- 前向验证"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 问题框架\n",
    "有许多方法可以利用和探索家庭用电数据集。\n",
    "\n",
    "在本教程中，我们将使用数据探索一个非常具体的问题。那是：\n",
    "\n",
    "“考虑到最近的功耗，未来一周的预期功耗是多少？”\n",
    "\n",
    "这就需要一个预测模型来预测接下来七天内每天的总有功功率。\n",
    "\n",
    "从技术上讲，给定多个预测步骤，对问题的这种构架称为多步时间序列预测问题。利用多个输入变量的模型可以称为多元多步时间序列预测模型。\n",
    "\n",
    "这种类型的模型可能有助于家庭规划支出。在供电方面，这对于规划特定家庭的电力需求也可能会有所帮助。\n",
    "\n",
    "数据集的这种构架还表明，将每分钟的功耗观察结果降采样为每日总计将很有用。鉴于我们对每天的总电量感兴趣，所以这不是必需的，但很有意义。\n",
    "\n",
    "我们可以使用pandas DataFrame上的resample（）函数轻松实现此目的。使用参数' D ' 调用此函数可将按日期时间索引的加载数据按天分组（请参阅所有偏移别名）。然后，我们可以计算每天的所有观测值的总和，并为八个变量中的每个变量创建一个新的每日耗电量数据集。\n",
    "\n",
    "下面列出了完整的示例。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1442, 8)\n",
      "            Global_active_power  Global_reactive_power    Voltage  \\\n",
      "datetime                                                            \n",
      "2006-12-16             1209.176                 34.922   93552.53   \n",
      "2006-12-17             3390.460                226.006  345725.32   \n",
      "2006-12-18             2203.826                161.792  347373.64   \n",
      "2006-12-19             1666.194                150.942  348479.01   \n",
      "2006-12-20             2225.748                160.998  348923.61   \n",
      "\n",
      "            Global_intensity  Sub_metering_1  Sub_metering_2  Sub_metering_3  \\\n",
      "datetime                                                                       \n",
      "2006-12-16            5180.8             0.0           546.0          4926.0   \n",
      "2006-12-17           14398.6          2033.0          4187.0         13341.0   \n",
      "2006-12-18            9247.2          1063.0          2621.0         14018.0   \n",
      "2006-12-19            7094.0           839.0          7602.0          6197.0   \n",
      "2006-12-20            9313.0             0.0          2648.0         14063.0   \n",
      "\n",
      "            sub_metering_4  \n",
      "datetime                    \n",
      "2006-12-16    14680.933319  \n",
      "2006-12-17    36946.666732  \n",
      "2006-12-18    19028.433281  \n",
      "2006-12-19    13131.900043  \n",
      "2006-12-20    20384.800011  \n"
     ]
    }
   ],
   "source": [
    "# resample minute data to total for each day\n",
    "from pandas import read_csv\n",
    "# load the new file\n",
    "dataset = read_csv('household_power_consumption.csv', header=0, infer_datetime_format=True, parse_dates=['datetime'], index_col=['datetime'])\n",
    "# resample data to daily\n",
    "daily_groups = dataset.resample('D')\n",
    "daily_data = daily_groups.sum()\n",
    "# summarize\n",
    "print(daily_data.shape)\n",
    "print(daily_data.head())\n",
    "# save\n",
    "daily_data.to_csv('household_power_consumption_days.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 评估指标\n",
    "预测将由七个值组成，每个星期的每一天都有一个值。\n",
    "\n",
    "对于多步预测问题，通常需要分别评估每个预测的时间步。这很有帮助，原因如下：\n",
    "\n",
    "- 在特定的交货时间（例如+1天与+3天）的性能进行评估。\n",
    "- 根据他们在不同交货时间的性能来对模型进行对比（例如，擅长+1天的模型与擅长+5天的模型）。\n",
    "\n",
    "总功率的单位是千瓦，因此具有相同单位的误差度量将很有用。均方根误差（RMSE）和均方根绝对误差（MAE）均符合此要求，尽管RMSE更常用，并将在本教程中采用。与MAE不同，RMSE更加惩罚预测错误。\n",
    "\n",
    "此问题的性能指标将是从第1天到第7天的每个交付周期的RMSE。\n",
    "\n",
    "作为捷径，使用单个分数总结模型的性能可能会很有用，以帮助模型选择。\n",
    "\n",
    "可以使用的一种可能得分是所有预测天的RMSE。\n",
    "\n",
    "下面的函数validate_forecasts（）将实现此行为，并基于多个7天的预测返回模型的性能。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# evaluate one or more weekly forecasts against expected values\n",
    "def evaluate_forecasts(actual, predicted):\n",
    "\tscores = list()\n",
    "\t# calculate an RMSE score for each day\n",
    "\tfor i in range(actual.shape[1]):\n",
    "\t\t# calculate mse\n",
    "\t\tmse = mean_squared_error(actual[:, i], predicted[:, i])\n",
    "\t\t# calculate rmse\n",
    "\t\trmse = sqrt(mse)\n",
    "\t\t# store\n",
    "\t\tscores.append(rmse)\n",
    "\t# calculate overall RMSE\n",
    "\ts = 0\n",
    "\tfor row in range(actual.shape[0]):\n",
    "\t\tfor col in range(actual.shape[1]):\n",
    "\t\t\ts += (actual[row, col] - predicted[row, col])**2\n",
    "\tscore = sqrt(s / (actual.shape[0] * actual.shape[1]))\n",
    "\treturn score, scores"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 训练与测试集\n",
    "我们将使用前三年的数据来训练预测模型，并使用最后的一年来评估模型。\n",
    "\n",
    "给定数据集中的数据将分为标准周。这些周从星期日开始，到星期六结束。\n",
    "\n",
    "这是使用模型的选定框架的一种现实且有用的方法，其中可以预测未来一周的功耗。这对于建模很有帮助，在建模中，可以使用模型来预测特定的一天（例如，星期三）或整个序列。\n",
    "\n",
    "我们将数据分为标准周，从测试数据集开始倒推。\n",
    "\n",
    "数据的最后一年是2010年，2010年的第一个星期日是1月3日。数据于2010年11月中旬结束，数据中最接近的最后一个星期六为11月20日。这给出了46周的测试数据。\n",
    "\n",
    "split_dataset（）将每日数据分为训练和测试集，并将每个数据组织为标准周。\n",
    "\n",
    "特定的行偏移量用于使用数据集的知识来拆分数据。然后，使用NumPy split（）函数将拆分后的数据集组织成每周数据。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(159, 7, 8)\n",
      "(3390.46, 1309.2679999999998)\n",
      "(46, 7, 8)\n",
      "(2083.4539999999984, 2197.006000000004)\n"
     ]
    }
   ],
   "source": [
    "# split into standard weeks\n",
    "from numpy import split\n",
    "from numpy import array\n",
    "from pandas import read_csv\n",
    " \n",
    "# split a univariate dataset into train/test sets\n",
    "def split_dataset(data):\n",
    "\t# split into standard weeks\n",
    "\ttrain, test = data[1:-328], data[-328:-6]\n",
    "\t# restructure into windows of weekly data\n",
    "\ttrain = array(split(train, len(train)/7))\n",
    "\ttest = array(split(test, len(test)/7))\n",
    "\treturn train, test\n",
    " \n",
    "# load the new file\n",
    "dataset = read_csv('household_power_consumption_days.csv', header=0, infer_datetime_format=True, parse_dates=['datetime'], index_col=['datetime'])\n",
    "train, test = split_dataset(dataset.values)\n",
    "# validate train data\n",
    "print(train.shape)\n",
    "print(train[0, 0, 0], train[-1, -1, 0])\n",
    "# validate test\n",
    "print(test.shape)\n",
    "print(test[0, 0, 0], test[-1, -1, 0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "训练数据集确实具有159周的数据，而测试数据集则具有46周的数据。\n",
    "\n",
    "我们可以看到，第一行和最后一行的训练和测试数据集的总有功功率与特定日期的数据匹配，这些特定日期我们定义为每组标准周的界限。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 前向验证\n",
    "将使用称为前向验证的方案评估模型。\n",
    "\n",
    "这是需要模型进行一周预测的地方，然后该模型可以使用该周的实际数据，以便可以将其用作在随后一周进行预测的基础。这对于在实践中如何使用模型是现实的，并且对模型有好处，使他们能够利用最佳的可用数据。\n",
    "\n",
    "我们可以在下面通过分离输入数据和输出/预测数据来说明这一点。\n",
    "\n",
    "Input, \t\t\t\t\t\tPredict   \n",
    "[Week1]\t\t\t\t\t\tWeek2   \n",
    "[Week1 + Week2]\t\t\t\tWeek3   \n",
    "[Week1 + Week2 + Week3]\t\tWeek4   \n",
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "下面的一个名为validate_model（）的函数提供了用于评估此数据集上的预测模型的前向验证方法。\n",
    "\n",
    "提供scikit-learn模型对象作为函数的参数，以及训练和测试数据集。提供了一个附加参数n_input，用于定义该模型将用作输入以进行预测的先前观测值的数量。\n",
    "\n",
    "稍后将介绍scikit学习模型如何拟合和进行预测的细节。\n",
    "\n",
    "然后使用先前定义的valuate_forecasts（）函数针对测试数据集评估模型做出的预测。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 递归多步预测\n",
    "大多数预测建模算法会将一些观测值作为输入并预测单个输出值。\n",
    "\n",
    "因此，它们不能直接用于进行多步时间序列预测。\n",
    "\n",
    "这适用于大多数线性，非线性和整体机器学习算法。\n",
    "\n",
    "可以将机器学习算法用于进行多步时间序列预测的一种方法是递归使用它们。\n",
    "\n",
    "这涉及对一个时间步长进行预测，采用该预测值，并将其作为输入输入模型中，以预测随后的时间步长。重复该过程，直到预测了所需的步骤数为止。\n",
    "\n",
    "- prediction(t+1) = model(obs(t-1), obs(t-2), ..., obs(t-n))   \n",
    "- prediction(t+2) = model(prediction(t+1), obs(t-1), ..., obs(t-n))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在本节中，我们将使用递归模型进行多步预测，开发一种测试工具，用于拟合和评估scikit-learn中提供的机器学习算法。\n",
    "\n",
    "第一步是将准备好的训练数据以窗口格式转换为单个单变量序列。\n",
    "\n",
    "下面的to_series（）函数会将每周多变量数据列表转换为每天总功耗的单个单变量序列。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接下来，需要将日常用电的顺序转换为适合于监督学习问题的输入和输出。\n",
    "\n",
    "该预测将是前几天消耗的总功率的某种函数。我们可以选择以前的天数作为输入，例如一两个星期。总会有一个输出：第二天消耗的总功率。\n",
    "\n",
    "该模型将适合先前时间步长的真实观测值。我们需要遍历每日耗电的顺序并将其分为输入和输出。这称为滑动窗口数据表示。\n",
    "\n",
    "下面的to_supervised（）函数实现了此行为。\n",
    "\n",
    "它以每周数据列表作为输入，并以先前的天数作为创建的每个样本的输入。\n",
    "\n",
    "第一步是将历史记录转换为单个数据系列。然后枚举该系列，每个时间步创建一对输入和输出对。考虑到前几天的观察结果，这一问题的框架将使模型能够学习预测一周中的任何一天。该函数返回准备训练模型的输入（X）和输出（y）。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Defined 10 models\n",
      "lars: [388.388] 411.0, 389.1, 338.0, 370.8, 408.5, 308.3, 471.1\n",
      "en: [469.337] 452.2, 451.9, 435.8, 485.7, 460.4, 405.8, 575.1\n",
      "ridge: [387.659] 407.9, 388.6, 337.5, 371.2, 407.0, 307.7, 471.7\n",
      "huber: [392.465] 412.1, 388.0, 337.9, 377.3, 405.6, 306.9, 492.5\n",
      "sgd: [384.637] 398.1, 389.1, 333.7, 371.1, 403.4, 307.0, 468.7\n",
      "pa: [391.397] 409.5, 384.8, 333.3, 378.8, 403.2, 303.5, 497.2\n",
      "lr: [388.388] 411.0, 389.1, 338.0, 370.8, 408.5, 308.3, 471.1\n",
      "llars: [388.406] 396.1, 387.8, 339.3, 377.8, 402.9, 310.3, 481.9\n",
      "ranscac: [440.001] 486.7, 436.6, 336.2, 408.1, 449.2, 381.1, 548.8\n",
      "lasso: [386.838] 403.6, 388.9, 337.3, 371.1, 406.1, 307.6, 471.6\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXd4lFX2xz/vTCaZ9N4TSCghIQkJCU0QRKOgIGCvWHCt2Na1rrqoqKvuj13WtaOsiKDCooCAhSIqIC2hBEiFEEjvPZnJlPf3x5tAQiakN7if55lnMvct90wIZ86ce+73SLIsIxAIBIILF1VfGyAQCASCnkU4eoFAILjAEY5eIBAILnCEoxcIBIILHOHoBQKB4AJHOHqBQCC4wBGOXiAQCC5whKMXCASCCxzh6AUCgeACx6qvDQDw8PCQg4KC+toMgUAgGFAkJCQUy7Ls2dZ5/cLRBwUFER8f39dmCAQCwYBCkqRT7TlPpG4EAoHgAkc4eoFAILjAEY5eIBAILnD6RY5eIAAwGAxkZ2ej0+n62pR+jVarJSAgAI1G09emCAYIwtEL+g3Z2dk4OjoSFBSEJEl9bU6/RJZlSkpKyM7OJjg4uK/NEQwQROpG0G/Q6XS4u7sLJ38eJEnC3d1dfOu5EEhcDYsj4FUX5TlxdY9NJSJ6Qb9COPm2Eb+jC4DE1bDhCTDUKa8rspTXAKNu6fbpREQvEAgEvc22hWedfCOGOmW8BxCOXiBogoODQ1+bILgYqMju2HgXEakbwYBl3cEc/u/nVHLL6/BzseXZ6SO4brR/r80vyzKyLKNSiXhJ0EGcA5R0jaXxHkD8hQoGJOsO5vDX746QU16HDOSU1/HX746w7mBOt9y/urqauLg4YmJiiIyMZP369QBkZmYSFhbG/PnziYmJISsri3vvvZeIiAgiIyNZvHhxt8wvuMCZ+kLLMY0txC3okelERC/ol7y24RhJuZWtHj94upx6k7nZWJ3BxHNrEvl632mL14z0c+KVWeHtml+r1bJ27VqcnJwoLi5mwoQJzJ49G4DU1FQ+//xzPvzwQxISEsjJyeHo0aMAlJeXt+v+goscqSHGtveEmmIlko9b0CMLsSAcvWCAcq6Tb2u8o8iyzIsvvsjvv/+OSqUiJyeHgoICAAYPHsyECRMAGDJkCBkZGTz++OPMnDmTadOmdcv8gguc/UvBIwQe3Qe9UEUlHL2gX9JW5D3p7V/IKa9rMe7vYsuqhy7p8vwrV66kqKiIhIQENBoNQUFBZ2rX7e3tz5zn6urK4cOH+fnnn/nggw9YvXo1//3vf7s8v+ACJvcQ5MTD1e/0ipMHkaMXDFCenT4CW4262ZitRs2z00d0y/0rKirw8vJCo9Gwfft2Tp2yrAZbXFyM2Wzmxhtv5PXXX+fAgQPdMr/gAiZ+KWjsIOq2XptSRPSCAUljdU1PVd3ceeedzJo1izFjxhAdHU1oaKjF83Jycpg3bx5ms5Iyeuutt7plfsEFSl05HFkDkTeBrUuvTSscvWDAct1o/24vp6yurgbAw8OD3bt3WzynceEVICoqSkTxgvZz+Bsw1MKYP/XqtCJ1IxAIBL2BLCtpG/9Y8Ivu1amFoxcIBILeIHMHFKfB2Pt7fWrh6AUCgaA32L8UtC4Qfn2vTy0cvUAgEPQ0VfmQshFGz1V2wPYywtELBAJBT3NgOZiNMOa+PpleOHqBQCDoSUxGSFgGQy4H96F9YkK7HL0kSZmSJB2RJOmQJEnxDWOvSpKU0zB2SJKkGU3O/6skScclSUqVJGl6TxkvEAgE/Z70n6Eyp08WYRvpSB395bIsF58ztliW5UVNByRJGgncBoQDfsBWSZJCZFk2dc1UgeAcElcrjRoqsntcFEog6DT7PwMnfwi5utlwXv56Mk4sQqfPQ2vjy5Chz+DrM6dHTOiJ1M0c4BtZlvWyLJ8EjgPjemAewcVMYyu2iixAPtuKrRv6bq5YsYJx48YRHR3NQw89hMlkwsHBgZdeeomoqCgmTJhwRuBMIDgvJSfgxC8Qcw+oz8bVefnrSUl5CZ0+F5DR6XNJSXmJvPz1PWJGeyN6GdgsSZIMfCLL8pKG8cckSbobiAeelmW5DPAH9jS5NrthTCBoPz++APlHWj+evR9M+uZjhjpY/xgkfGH5Gp9IuObt806bnJzMqlWr2LVrFxqNhvnz57Ny5UpqamqYMGECb775Js899xyffvopL7/8cgfflOCiI+FzUFlBzN3NhjNOLMJsbi7KZzbXkXFiUY9E9e119JNkWc6VJMkL2CJJUgrwEfA6yofA68A/gfsAS3Js8rkDkiQ9CDwIMGjQoE6YLrioOdfJtzXeTrZt20ZCQgJjx44FoK6uDi8vL6ytrbn22msBiI2NZcuWLV2aR3ARYKiDgysgdCY4+TY7pNPnWbyktfGu0i5HL8tybsNzoSRJa4Fxsiz/3nhckqRPgY0NL7OBwCaXBwC5Fu65BFgCMGbMmBYfBIKLnDYibxZHtNKKLRDmber0tLIsc88997QQJ1u0aBFSg6SsWq3GaDR2eg7BRcKxdVBXZlHXRmvj25C2aTneE7SZo5ckyV6SJMfGn4FpwFFJkppadD3QqPT0PXCbJEk2kiQFA8OBfd1rtuCiJ25By40n3dCKLS4ujjVr1lBYWAhAaWlpqxLFAsF52f8ZuA+H4CktDgUP+UuLMZXKliFDn+kRU9oT0XsDaxuiGSvgK1mWf5Ik6UtJkqJR0jKZwEMAsiwfkyRpNZAEGIFHRcWNoNtprK7p5qqbkSNH8sYbbzBt2jTMZjMajYYPPvigGwwWXFTkHW5oLvK2xeYiGisH5VnjhsFQ1uNVN5Is933WZMyYMXJ8fHxfmyHoY5KTkwkLC+trMwYE4nfVz/m+oQLs6RSLuvMHD95DTe1xJl7yGypV59XiJUlKkGV5TFvniZ2xAoFA0J3oKuDI/1ptLlJbe5LSsp34+93eJSffEYSjFwgEgu6ksbnIWMvNRbJzvkKSNPj53dprJglHLxAIBN2FLCtyxP6x4De6xWGTqY68vDV4eU7Hxsaz18wSjl4gEAi6i8ydUJzaaqvAgoINGI2VBATc1atmCUcvEAgE3UV8Q3ORiBtaHJJlmezsFTjYj8DZObZXzRKOXiAQCLqDqnxI3tBqc5HKykNUVR/DP2Dumc13vYVw9ALBeZgxYwbl5eUtxl999VUWLVpk4QrBRcuBL8/bXCQ7ZwVqtQM+3j1TK38+eqe2RyDoATZlbOLdA++SX5OPj70PT8Y8ycwhM7vt/rIss3HjRlQqEQ8J2uBMc5GpFpuL1NeXUFDwA/7+t2JlZd/b1omIXjAw2ZSxiVf/eJW8mjxkZPJq8nj1j1fZlNF5nRuAzMxMwsLCmD9/PjExMajVaoqLlTYMb775JiNGjODKK68kNTX1zDX79+9n1KhRXHLJJTz77LNEREQAYDKZePbZZxk7diyjRo3ik08+6ZJtgn5M+s9Qmd1qc5HcvDXIcj0B/nN72TAFEdEL+iXv7HuHlNKUVo8nFiVSb65vNqYz6ViwawFr0tZYvCbULZTnxz3f5typqal8/vnnfPjhhwQFBQGQkJDAN998w8GDBzEajcTExBAbqyyozZs3jyVLljBx4kReeOGFM/dZunQpzs7O7N+/H71ez6RJk5g2bRrBwcFt2iAYYOxfCo5+EHJNi0OybCInZyWuLhOwtx/WB8aJiF4wQDnXybc13hEGDx7MhAkTmo3t2LGD66+/Hjs7O5ycnJg9ezYA5eXlVFVVMXHiRADuuOOOM9ds3ryZ5cuXEx0dzfjx4ykpKSE9Pb3L9gn6GaUZcGIbxDZvLtJISclv6HQ5+Af0TTQPIqIX9FPairynrZlGXk1L7W5fe18+v/rzLs1tb285h2qpUuJ8WlGyLPPee+8xfbpom3xBE/9fkNRKFykLZOeswMbaG0+PK3vZsLOIiF4wIHky5km0am2zMa1ay5MxT/bIfFOmTGHt2rXU1dVRVVXFhg0bAHB1dcXR0ZE9e5Smat98882Za6ZPn85HH32EwWAAIC0tjZqamh6xT9BHGHRwcKXF5iIAtbWnKCn5HT//21CpNH1goIKI6AUDksbqmp6sumlKTEwMt956K9HR0QwePJjJkyefObZ06VIeeOAB7O3tmTp1Ks7OzgDcf//9ZGZmEhMTgyzLeHp6sm7duh6xT9BHJK2DutJWF2Fzcr9CktT496KujSWETLGg3zBQpXerq6txcFD0xd9++23y8vJ49913e3TOgfq7uuD47EqoK4fH9rfQnTeZdOzcNQk3t4lERrzXI9O3V6ZYRPQCQRfZtGkTb731FkajkcGDB7Ns2bK+NknQG+QdVprUT3/LYnORgsKNGI3lBPjf2QfGNUc4eoGgi9x6663cemvffjUX9AH7l4KVLUTfbvFwTvZK7O2H4+IyvpcNa4lYjBUIBIKOcqa5yI1g69ricGVlIpVVifj739nrujaWEI5eIBAIOsrhVUpzkVbkiLOzV6BW2+Prc10vG2YZ4egFAoGgI8gy7P8M/GLAP6bFYYOhjILCjfj4XIeVlWMfGNgS4egFAoGgI5zapTQXaaVVYG7eGsxmfb9YhG1EOHqBoAmZmZlnRMnaw7Jly3jsscd60CJBv2P/UtA6Q7il5iJmcrK/wsVlHA4OI/rAOMsIRy8YsFRs2ED6FXEkh40k/Yo4Khp2qw4kTCZTX5sg6AhVBZD8PUTPBWu7FodLSn+nTne6X0XzIBy9YIBSsWEDeX9bgDE3F2QZY24ueX9b0C3O3mQy8cADDxAeHs60adOoq6tj6tSpNG7qKy4uPqNqCZCVlcXVV1/NiBEjeO21186Mr1ixgnHjxhEdHc1DDz10xqk7ODiwYMECxo8fz+7du7tsr6AXObj8vM1FcrJXYm3tiafntF427PyIOnpBvyT/739Hn9y6THHd4cPI9c2VKmWdjryXXqZ89f8sXmMTForPiy+2OXd6ejpff/01n376Kbfccgvffvvtec/ft28fR48exc7OjrFjxzJz5kzs7e1ZtWoVu3btQqPRMH/+fFauXMndd99NTU0NERERLFy4sE1bBP0Iswnil0HwZeDRUm64ri6L4pLtBAU9ikpl3fv2nQfh6AUDknOdfFvjHSE4OJjo6GgAYmNjyczMPO/5V111Fe7u7gDccMMN7Ny5EysrKxISEhg7diwAdXV1eHl5AaBWq7nxxhu7bKegl0lraC5y9VsWD+fkfI0kqfD3u62XDWsb4egF/ZK2Iu/0K+KUtM05WPn5MfjL5V2a28bG5szParWauro6rKysMJvNAOh0umbnn7shRpIkZFnmnnvu4a23WjoFrVaLWq3uko2CPiB+KTj6wogZLQ6ZTHpy81bj4XElWm1LFcu+pl05ekmSMiVJOiJJ0iFJkuIbxtwkSdoiSVJ6w7Nrw7gkSdJ/JEk6LklSoiRJLQtNBYIu4vXUn5G0zWWKJa0Wr6f+3CPzBQUFkZCQAMCaNc07WG3ZsoXS0lLq6upYt24dkyZNIi4ujjVr1lBYWAhAaWkpp06d6hHbBL1AaQYc3wqx91psLlJY+AMGQ1m/W4RtpCOLsZfLshzdRCntBWCbLMvDgW0NrwGuAYY3PB4EPuouYwWCRpxnzcL39YVY+fmBJGHl54fv6wtxnjWrR+Z75pln+Oijj5g4ceKZHrKNXHrppdx1111ER0dz4403MmbMGEaOHMkbb7zBtGnTGDVqFFdddRV5eS0bpQgGCPGfNzQXudvi4eycldjZDcHVdWIvG9Y+2iVTLElSJjBGluXiJmOpwFRZlvMkSfIFfpVleYQkSZ80/Pz1uee1dn8hUywAIb3bEcTvqhcx6OBfYRA0CW5d0eJwZdVR9u+fQ8jwvxEYeG+vmtbdMsUysFmSJBn4RJblJYB3o/NucPZeDef6A1lNrs1uGBPhjEAgGHi01VwkeyUqlS0+Pi03UJ2P5B3b2fHNcqpKinF092DybXcTNvny7rC4Be119JNkWc5tcOZbJElqve4NLEm1tfjaIEnSgyipHQYNGtROMwQCgaCX2b8U3IcpZZXnYDBUkF/wPT4+16HROLX7lsk7trN5yfsY6/UAVBUXsXnJ+wA94uzblaOXZTm34bkQWAuMAwoaUjY0PBc2nJ4NBDa5PABoUR4hy/ISWZbHyLI8xtPTs/PvQCAQCHqKvETI3qdskLIgN5yX9y1ms44A/7kduu2Ob5afcfKNGOv17PimaxVjrdGmo5ckyV6SJMfGn4FpwFHge6Cx7fk9wPqGn78H7m6ovpkAVJwvPy8QCAT9lvjG5iJ3tDgky2ayc1bi7ByLo2PH1kuqSoo7NN5V2pO68QbWNtQKWwFfybL8kyRJ+4HVkiT9CTgN3Nxw/g/ADOA4UAvM63arBQKBoKfRVUDi/yDCcnOR0tJd1NVlMiT4yQ7f2sHOnuqaaovjPUGbjl6W5QwgysJ4CRBnYVwGHu0W6wQCgaCvOLwKDDUw1rKuTXbOCjQaN7y8pnf41sPySjnkqGmWDlKZzYTklXba3PMhRM0Egm4iKCioRY29YIAiy0raxm80+Me2OKzT5VJc/Av+freiUtlYuMH5UZeWgiRhbTCCLKOtNxCZVYRPZnZ3WN8CIYEgGLCk7c1n9/oTVJfqcXCz4ZI5QwkZ79PXZgkuBE79AUUpMPt9i4dzcr4CwN+/Ze6+PeT5eGBtNHFF0qlm0baVn1+n7tcWA9bR5+WvJ+PEInT6PLQ2vgwZ+gy+PnP62ixBL5G2N5/tK1Mw1iv6M9WleravVKp+u+Lsa2pquOWWW8jOzsZkqOdvT/4JRzsb/rJwMR6eXsSMGUdGRgYbN26kpKSE22+/naKiIsaNG0d7Nh8KBgj7P1Oai0S0FJ8zm/Xk5K7Gw+MKtNqOO2Z9bQ0FdjYEFJY2c/I9KeExIB19Xv56UlJewmyuA0CnzyUl5SUA4ewvEHasTqM4q+ViVSMFJyswGZs7VmO9mV++TObYzpZiZwAegQ5MviXkvPP+9NNP+Pn5sel/X0JFFhUVFURccQu/f/cZwYMDuf2Js9LCr732GpdeeikLFixg06ZNLFmypAPvUNBvqS6E5A0w7gGLzUUKC3/GYCjpcEllI2l7d2Eym/Avq0Ll7Iy5shIrX1+8nvpzj0l4DMgcfcaJRWecfCNmcx0ZJxb1kUWC3uZcJ9/WeHuJjIxk69atPP/cs+zYE8/J07kMGexP8CB/kM3cPuvsZpbff/+duXOV/+wzZ87E1bVlZYZgAHJgOZgNrTYXyc5Zga3tYNzcJnXq9sk7fsXR2gYXEwzbuoWw5CSG/7Ktx5w8DNCIXqe3XJbf2rhg4NFW5P3Fi7uoLtW3GHdws+H6pzsvmBoSEkJCQgI/rPyQv771PldNGd/8BLOx2ctzJYoFAxyzCRKWNTQXGd7icFVVMhUVCQwf9hKS1PE4ubK4iKykI4SUVlN4zdX89NlnVFRU4OzsTFxcHKNGjeqGN9GSARnRa21a03uWyc5eKXKlFwGXzBmKlXXzP18raxWXzBnapfvm5uRgJ9cy98aZPPPwXfwRf5iMUzlkZinpoFXfb4b6WqgrY8rkyaxcuRKAH3/8kbKysi7NLegHpG+GiiwY+yeLh7NzVqBSafH17VzjmJRdv4EsY7LSssPOjoqKCgAqKirYsGEDiYmJnTb9fAzIiH7I0Gea5egBVCottraDSU1bQEnp74SFvoW1tVsfWinoSRoXXLu16sZs5MjurTy74A1UKjUaKzUfvfUieQVFXH3nY3i4uTAudjRQAGWZvPLwzdz++MvEfLuGyy6bKjSbLgT2LwUHH4vNRYzGKvLz1+PjPRuNxrlTt0/e+Stuag1poyIwNjSyacRgMLBt27YeieoHpKNvXHA9t+rGx3sWWdlfcPz4P9i7bybhIxd1Oo8m6P+EjPfpvnJKfRWUnWL6xCim790JDl5QVwZVeVRXlpOyayOygw+PPvs3xlwaB+7DcLctZfPK90A2g9qaxQufA1uH7rFH0PuUnlSai1z2PKg1LQ4rujZ1+Ad0rrlI0amTFJ/OJDynmFMhLfagApyJ8LubAenoQXH2lipsBgXOw9VlHEePPcXBQ/cweNCDDBny537XrFfQT5DNUJUP1QWgtgGPkLOVFnZuYOfGp18t5osvvqC+vp7Ro0fz0MMPg40d2DiCOUDZKl9XqtyjugA0tmDrpmybt+AwBP2UhM9BUkHsPS0OybJMds5KnJyicXKM6NTtk3ZsR5IkfMuqcLK3p7KmpsU5zs6d+6bQFgMyR98Wjo7hjBu7Dj+/Wzl1+hPiE26mtvZkX5sl6G8YdVCcrjhnO3fwHGGxnO6pp57i0KFDJCUlsXLlSuzsmpyjUisfCO7DwDsCnPwVUe7KHCg4CiUnoLZUWeQT9F8MOjjwJYy4Bpxa1saXlf1BbW1Gp0sqzWYTyTt/xdsg4xwVzZTLW0oRazQa4uJaqMp0CxekowdQq+0IC32TyIgPqavLYt/+2eTmrhELtQJle3ttCRSlglEPrsHgMkhx2l1BrVFSPl6h4Bmq/Gyog/JTitMvO6WkiMTfYP8jaf15m4tk56xs0LVpmbtvD1nHjlBTVopvdj7O11+HyaR88Ds4KKk+Z2dnZs2a1WNVNwM2ddNevLym4+Q0imNJT5Oc8jwlpb8ROuLNDjUJEFxAmI1QngW6crB2AJfBYNUDaT2NLWj8wdEP6quVqF5XrjgTlUb5FmDrqpwn6Hvil4LbUIvNRXS6PIqLtzIo8H7U6o7r2oBSO2+lUuFdb8Jx+nTiv/wSf39/Hnjgga5a3i4ueEcPoNX6EjP6S06d+pSMk4uprDhEePhiXFzabLUouJDQVyvRtckAjr7g4G2xmUS3IklKLr9ZPr9M5PP7E/lHIGsvTHsTVC2THDm53yDLZvz9b+/U7Q16Hel7d+FbUYPLlVeRU1ZGUVER+cNG8Z8f4inWqvHQmbiyroTFN3ZcCbM9XLCpm3ORJDVBQQ8TG7saSWVFwoHbycj4N+ZzNsAILkBkM1TmQkm68tpjODj69LyTP5cz+fyhZ/P50CSff1zk8/uC/UvBSmuxuYjZXE9u7irc3adiaxto4eK2OZGwj3pdHb6FZbhcfx3x8fEc9w7kR59gimytkCWJIlsrvnP24qlvf+7qu7HIRePoG3F2imLc2A34+MzhZOZ7HDh4O3V1WW1fKBiYGPVnF1xt3ZTcuXXPNHfoEI35fM/GfL63siAo8vm9i64SElcr4mV2LffdFBVtpr6+iICAzi3CgtIf1lZS4e3oghwRQVJSEvGDw9GrmwcaerXEVlv3Ts9zPi6K1M25WFk5ED5yEe5uU0hJ/Rt7911L6IjX8fGZ3demCTpA8o7t7PhmOVUlxTi6ezD5trubN1auLVV2OSKBa5DFLkH9Ao2t8nD0VfL5daVQ1zSf76p8SIl8fveT2NhcpLWdsCux1Q7C3W1Kp25fW1lB5qEDBBWW4jLnJg4fOYLJZKJCaznXX6ztYkFAK1x0EX1TfHxmM37cRhwcQjiW9BTHkp7BaGxdMVHQf0jesZ3NS96nqrgIZJmq4iI2L3mf5B3blQXXskwlOtbYKmWT7XTymZmZhIaGcs899zBq1ChuuukmamtrWbhwIWPHjiUiIoIHH3ywZ6q3GvP5LoPBO1L5cNLYKmqKRSlQmKL8bDJ0/9wXI7KspG18oy02F6muTqW8fB/+AXd0StcGIHX3DsxmE36lVThfp6RtBg8ejKfObPF8r1bGu8pFGdE3xdY2kJjRX5OZ+QEnM9+noiKB8PB/4+xkeeeaoHfYvmwJhacyWj2el5aKydjc4Rnr9fz88bsk/rT2zG5V1NbQ8A3Za/AQLr/3wTbnTk1NZenSpUyaNIn77ruPDz/8kMcee4wFCxYAcNddd7Fx40Zm9aDaICqV8uFk66o49roy5VGZozyqqyDxCITO7B+pqIHI6d1QlAyz37N4ODtnJSqVNX6+N3V6iuQd23EyyXiHR5JtMlFWVsYVV1xB7K5qfhzR/Bua1igz40id0nG7m7moI/pGVCorhgx5ktiYr5HNRhISbiEz8yNkWSyK9VfOdfJnxxsW1zV2StlkJ9ZbAwMDmTRJkc6YO3cuO3fuZPv27YwfP57IyEh++eUXjh071lnTO86ZfP4I8AxT8vlmA3z3ACwKgbUPw4ntYhG3o+z/DGycIaKlI1d0bdbh7XUtGk3nUn5leTnkpafiV1B6ZhHWzs6OsLAw8l1U2BvMeNWZkWQZ71oztx+oIbSiZwoELvqIvikuLmMYN24TKakvcSJjESWlOwgf+U+02tbUMgU9RVuR95JH5ylpm3NwdHPl1tcXd2nz07nSw5IkMX/+fOLj4wkMDOTVV19Fp9N1+v5dQqMFjR84lsO9m+DwN8pmn8NfKzn+yJsh6jbwDu8b+wYK1YWQ9L2yQcrCbui8/HWYTDUEBNzV6SmSd/4KgF9dPVx6KSmffMLEiRPZnpXOIS9rbknTE3Ko9sz5VtYqLrkztNPznQ8R0Z+DRuNERPh/CAt9h6qqI+zdN5PCwp4peRJ0nsm33Y2VdfMFLStraybfcV+Xd7iePn2a3bt3A/D1119z6aWXAuDh4UF1dTVr1qzp0v27BUmCoEthzvvwTBrc9Dn4RsGeD+GjifDRpfDHe4qOT+JqWBwBr7ooz4mr+9r6vuc8zUVkWSYnZyWOjpE4OXVup6osyyT9vh33unq8rojjcGoqsiwTGxvLBwdOopIh9IQOB1flb9jBzYbL7wztsZ7HIqK3gCRJ+PndhItLLEePPcWRo/Px87uNkOEvoVa3/PQX9D5hk6ZATTE71q2jqqwCRw8PJt92T/Oqm87eOyyML774goceeojhw4fzyCOPUFZWRmRkJEFBQYwdO7Yb3kE3orGFiBuUR00xHP1WifQ3v6w8JJWyZgFKFdKGJ5SfR93Sdzb3JWeai0wBz5YNbsrL91FTk05Y2DudniIvPZWKwnwii8pxen4OCXv3MnToUOptVBxx9GZyvoFR4/yYdPP5G+x0F8LRnwc7u2DGxK4mI+PfnDq9hPLy/USE/xtHx5F9bdrFTX0NlGUSNmoYYZMWKfrh3bj5SaVS8fHHHzcbe+ONN3jjjTe6bY4wstQ1AAAgAElEQVQew94Dxj+kPIrS4LM40Fc2P8dQB9sWXryOPn2L8oE3zfK/Z3bOCqysnPH2urbTUyTv3I4KCLB1IMfNjcrKSq655hoWbN9Fjb0/Y9NriHyscxuwOoNI3bSBSmXNsGHPMTr6C4zGKvbH38jp0/9FlnumDEpwHmRZSUUUpymv3YcreWnRzs8yniHKpitLVGT3ri39if2fKcFB6MwWh/T6AoqKNuPnexNqtbZTtzcZjaTs/A3v8mo85swh4cABHB0dGTpsKPG4EVJhYoSXE04evbcvQjj6duLmNonx4zbh7j6F9ONvcujwfej1LRcDBT2EsV6RMKjKU0oOPUeATfc3+QgKCuLo0aPdft8+wznA8riVjbKh7GKjLFNpLhJ7j0VtoZzcVciyEX//zjUXAcg8fABdTTV+ZVVw5ZWkp6cTExPD+/F7yLW3ZWqGntFXDu7Cm+g4wtF3AGtrN0ZFfsyIkIWUl+9j774ZFBdv72uzLnzqypQNQ4Y6ZTORaxCoRNaxXcQtaLmjVq1RavOXXAa5h/rGrr4ivqG5SEzL5iJms4HcnG9wd5uCnV3nHXHyju1Ym2UGDQ/lSH4+kiQRExPDprx6nOvNDK8Cv+EuXXkXHWbgOvo+qiSQJImAgDsZO2YdNjZeHE68n9S0hZhM+s7fVFRFWMZsUjRfyjKVCNQz1KIeieA8jLoFZv0HnAMBSXme8yH8aYvy+/3vdDi4sq+t7B2MejjY0FzE2b/F4aLirejrC7pUUqmvreX4/t34llbifP31HDhwgJCQEBIrikh2cSfuVD2xVwa1KOHtadodFkmSpAbigRxZlq+VJGkZcBnQ2OTwXlmWD0nKO3gXZX9XbcP4gW61OnG1UjlgaGgOXpEF3z+hfL2PvBGQGvK2rTx3wy/ZwSGEMbHfceLEP8jKXkZ52R7Cw/+Ng0MHV9EtvZeLvSoCziy4YqpX8qmO3kokJug4o26x/Lf00O+wZh6snw858XD128oH6oVK0nql4UwrujY52SvQav1xd2+pSd9e0vfuwmQ04l9bT27wEGqOHiU2Npa/7k8FVz8iM/WE3Ovd6ft3lo58/30SSAaadux4Vpblc4uKrwGGNzzGAx81PHcf2xaedYyNGOvg+0eVR7to48PgvM/K9WpJIgQJNxcVSUGp7N9zDcOzJfyLVQ3aGK1fe+Z1RZaizdIUQx1sfeXidPSyrChNVuUrKQb34T2Si7+YyMtfT8aJRej0eWhtfBky9Bml37K9B8xdC78shF3vQl4i3LLcYrR7QbD/M3AbAsFTWxyqqTlOWfkehg55FiWm7RxJv2/DzmBk0OTL2Hr0CM7Ozrj6e3HkZBWXFBiJHuOHlXXPCJedj3Y5ekmSAoCZwJvAX9o4fQ6wXFZUn/ZIkuQiSZKvLMt5XTO1CeerGIh7BZAb5F1lpX9ns9dyQ03xOWOtPtPmeR6ymfGGOpKkQ6QOKqLEz50wXSTWssby/E2fy1rpZVuZC/8ZDYETYNB4GHSJ4vQsNEa4YDAbFU32+mrQuoBLYK/n4h0cHKiuvnCE7fLy15OS8hJmsxIY6fS5pKS8BKA4e7UVXLVQEfVaNx8+mQI3f67UmF9I5B89b3OR7JwVSJI1fn43d3qKqpJispKOMqykEvOj0zm5YwdXXHEFr/36B5V2foxLryXy8d4rqWxKe/8X/Rt4DnA8Z/xNSZIWANuAF2RZ1gP+QFOB9+yGsWaOXpKkB4EHAQYNGtQxq50DGuRnzx0PhMltfQ71DDZAtGwmK/sLjh//B3vtjhA+chFubpPOf+HpPZbfi9ZZ0TVJ/xkOf6WM2bpC4HjlMWgC+MUoW+IvBI6tgyp7cPNS+rfaurWZYqs5WEjlz5mYyvWoXWxwmh6E/WivbjfNZDKhVvd+FNYdZJxYdMbJN2I215FxYpHi6BsZOUf5e1t1JyyfA1e+ChOfuHBKV+Nbby5iNNaQl7cWb68ZWFt3Xg8+ZddvAAyytiXJYEClUhEdHc3Tvx0luMpEpKcTDq598/+1zfBQkqRrgUJZlhPOOfRXIBQYC7gBzzdeYuE2LTRdZVleIsvyGFmWx3h6enbMakuVBBpbZbwPkSQVgwLnMXbMGqysHDl46B6OH38Hs7m+9Ytaey8zFsHtX8GzJ+CxeJj9vlL3W3ICtr0Gn18DbwfCZ1cpux+TNyq7Igca+mpY/xj87x5Fe91zBNi5t8vJl3+XjqlcWQQ3lesp/y6dmoOF3WLWr7/+yuWXX84dd9xBZGRkt9yzt5FlGZ0+1+Ixnd7CF2zPEHjgFwibBVsWwOq7W6/DH0joKuHwKgi/weJifn7Bekym6i41FwFI2r4Fl1odPtfO5tDhw4SGhrIi+QhZDnZKSeVVvVtS2ZT2RPSTgNmSJM0AtICTJEkrZFlu/K3oJUn6HHim4XU20PT7SQBg+a+tszTmrrctVNI4zgGKw+wnOW1Hx3DGjV1PevqbnDq9hNKyP4gI/zd2dsEtT27rvUiS0vrOYzjENFQD1BQrX0NP71Ge936i6JoAuA87m+4JnKBc11+jstyDsOZPUJoBl/5FUWi0UiKe8g0nqM+tafXS+tOVYGoeP8gGM2Vr0qjZl2/xGms/e1xmDW23efv27ePo0aMEB1v4d+vnVNekk5a2sNXj1tatBFc2jnDzF8rf09ZX4NMUuHWlRamAAcOZ5iL3tzgkyzI52StwdAjHySm601MUnTpJcW42I0uryIuMpO733xkzZgwPH8zAwdGOkVUSPkOdu/IuukSbjl6W5b+iRO9IkjQVeEaW5bmNefeGKpvrgMZdJt8Dj0mS9A3KImxFt+bnG2mtkqCfoFbbEhr6Bm7uk0lOfpF9+2cTMvwVfH1vbFla1dH3Yu+hRPeNO/sMOsg7pOhrn94LqT/AoRXKMVs3Jc1zJt0zuu8rK8xm+OM/8MsbinO/ZwMET4bk5Pbfw9TiS+L5xzvBuHHjBpyTNxqryDj5H7Kzl6NW2+HjcwOFhT9gNuvOOa+G6upUHBxGtLyJJMGkJ8AvGv43Dz69HK77UEnvDDRkGeL/qwi++ce0OFxRkUB1TSqhoX/vUslj0o5fkWSZoUHD2HHyJO7u7hRp4JiLOzMy6xl7Ve+XVDalKytdKyVJ8kRJ1RwCHm4Y/wGltPI4SnnlvC5ZOMDx8pyOk+MojiU9TXLK85SU/kboiDfQaLrx012jVZz4oAnKa1lW+qSe3n028k/9QTmmtlZy+40Rf+B4sO+ZPpUWqcyFtQ/Byd8hbDbMetfi1+m2Iu+8t/edSds0Re1ig9dDnVMcPBd7+4HT0EOWzeTnr+X4iX9QX1+Cn98tDB3yNNbW7ri5Xdqs6sY/4G6ysz7nwME7GR29vHXtpuApSgnm6ruVx8QnlGIH9QDarHZ6NxQmKXsJLDja7OwvsbJyxMe7821EZbOZ5F+34FFVi2ruHWSlpzNt2jT+tS8Zs6sfUZl6ht/d+yWVTenQv5gsy78Cvzb8fEUr58hAe2scLwq0Wl9iRn/JqVOfknFyMZUVhxgZ/i9cXXpIBVGSlK/aniHKVm+A6qIGp9/g/Hd/qJTUgVLN0/hBETgB3If2TLoneSN8/5iycWX2ezD6rk7P4zQ9iPLv0pENZzWHJI0Kp+lB3WTswKGy8gipaa9RWXkQJ6fRRI36tJm8rq/PnOYLr4CX51UcODiXAwfnMjr6C5ycWlmHcPaHeT/AT39VvoXlHlQkkR06uK7WV+xfqjQXiWzZXESvL6Kw6GcCAuaiVndedyYr6Sg1VZWEVOtJsXdArVYTHBpCYmUmY4uMjB3jj1rTt9VyA+ijeWAjSWqCgh7G1e0Sjh37MwcO3EFw0KMEBT2GqjdKCB08Iexa5QFKrX7uQSXaP70HkjcouwYB7Dyap3t8o7qW7qmvgZ9fVKRhfaPhxqXgMaxLb6exuqY3qm76K/X1pZzIWERu7mqsrd0ZGfYPfHyub1d/Uzu7IGJjvubAwbkcPHQX0VHLcHZuJUdtZQPX/gsCxsDGpxTphFuWK6/7M9VFyiapsX+y2G4xN3cVsmwgoAu6NgBJv23Fymxm8IRLWZV0jPDwcP6xO54yOx8mHK8i4rFW9IZ6EeHoexlnpyjGjd1AWtprnMx8j9LSnXh5zSAr6/OWG1p6Eo0tDJ6oPEDJmxenQdaes84/ZaNyTG2j5DcbI/7Ace2XIsg7rCy4lhyHSU/C5S8rLf66AfvRXt3u2Btr6KdOncrUqVO79d7dhSzLGI3V7N5zByZTLYGB8xgS/ARWVudWP58fW9vABmd/JwcP3UN01FJcXM7jvKPvUDpXrZqrVH1d8w7Ezuu/i/0HG5uLtNwJazYbycn9GjfXSy0XSbQTQ72etN078S6vpnDcOOqPHWXMmDH838HTBNSYGO3pjL1z3+82Fo6+D7CycmDkyP/DzW0ySckvUFF58MyxFhtaeguVCrxClUfsvcpYVYGS5mlM+fzxHpgXK8c8RpzdyBU4XtlxKEmKpENjBZHWCXRVinzB3ethSOe3lgsUjMYadPpcjMYKHB0jCAlZgIP98E7fT6v1Iybmaw4enMuhw/OIGrUUV9dxrV/gGwUP/qb0q934FGTHw8x/tiwR7mvMJohfBkGTLVYMlZT8gl6fz4iQV7o0TUbCPgyGegarbEgsLcHLy4vNBTmcdHLgriO1xMzq/L9NdyIcfR/i4zOb48ffRl9f0Gzc4oaWvsDRG0bOVh4A9bWQe+BsWWfSeqUlG4C9Fzj5QcExJYoC0FWApIbLnhNOvouYzQb0+jwMhgpUKg0ajRvh4cu7pZJDa+NDzOivOHDwrgZnv+T8G/3s3OCO1fDbO8oj/wjc+qWiKtpfSN8CFadhmuUS0+zsFdjY+OLubnGpsd0c2/oTNgYjDtOuIS8vjxkzZvDqqQpsnbVEVkt4Bzu1fZNe4ALeTz8w0Ndb3uBjcUNLX2Ntp/QpnfIM3Pk/eC4THtkN1y6GoVco/+EbnXwjsgl2/KtPzL0QkGUzen0hNTVpGIyVWNt4YW8fglpt263lejY2XsTGrMTOdjCHEx+gpOT381+gUsPlL8LtqxSF0SVTFZ33/kL8UnDwhtCWXaJqajIoLduFv//tXVofq62sIPNYIn5l1Rz39UGj0YCXG0ddPJiaVc+Eq/pPaa5w9H2M1sa3Q+P9CpUKvEcqDZZv+ORsX9JzuZi7GXUBg7GKmpp09PoC1GoHHOyHo7Xxbtdia2ewtvZg9OgV2NkN5XDiQ+3rtTDianhwOzj5w4qb4Lf/U9Z7+pKyTCWij2mluUjOSiRJg5/frV2aJnX3DmRZJjAgiKSTJ4mMjOTd+FTq1RLRp/UM7UeFAcLR9zFDhj6DSnVuflPFkCF9o9nTJVrrZtTauMAiJpOe2tpM6mozAaVCxs5uMCpVzy/qWVu7ETP6SxwcRpB45BGKira0fZH7UEXfPvJm2P4GfHMH1JX3uK2tEv+5sl7UWFrcBJOplrz8b/Hyuhoba48uTXPsp0041OmpuXQKBoOBkVGRHLb3ZnSxkYljAlBb9R/32n8suUjx9ZlDaOibaG38AAkrKxfATH39AGxT2E81iAYKsmxGp8+npjYdk6kGGxsf7O2Hd7iapqtoNC4NG6nCOXL0MQoKf2z7Ims7uGEJXPN/cHyLspu24FjPG3suZ5qLzLAYYOQXbMBorCLAv2u6NuX5eRTkZuFfo+eYoR4/Pz8+PZZGsdaKS47rCJ/Sv4Ib4ej7Ab4+c5g0aQdxVxxnyuR4PD2ncSJjMVXVKX1tWsew1M1o1n/6tVTFuTg4KNr3mZmZRERE9MqcsixjMFRQXZNGvb4IjZUz9vYh2Nh49liapi00GidGRy/DySmKY8eeJD//+7YvkiQY/yDcu0lZuP80DhL/1/PGNiXpe6W5yJj7WhySZZns7BU4OITi7BzbtWl+3QKyjGvsBIqKixkzZgx/6O3wqTUzztMJO6fuKSHuLkTVTT9DkiRCR7zB3ooZJB37C2PHru2Vr+zdRi9qECUmJrJt2zYqKipwdnYmLi6OUaO6R/6go3RWythk0qHT52Iy1qBSa7G1C8TKqn9IL1hZORId9TmHEx/gWNLTyLIRX98b2r5w0ARFOuF/98J39yvdq6a9YTFf3u3s/wxcg2HI5S0OVVYepLo6iREjXu/SQrYsyxzb+jPu1Tqyhw3DpqSERGM96S6O3H6sjtjZ/aOksikiou+HWFu7Exb6NtU1qZzIWNzX5vRLEhMT2bBhAxUVSifLiooKNmzYQGJiYrfPlZmZyeTJk4mJiSEmJoY//vgDaCllXFNTw8yZM4mKiiIiIoJVq1a1ek9ZNqHT5VJTk47ZpEOr9cPebli/cfKNWFnZEx21FDfXS0hKfo7c3HZG6I7ecM/3MGE+7P0Yll2rdA3rSQqOKRv+xv7JcnOR7BWo1Q74eHetbDn/eBqVVRX4WGlJKyggKiqK7zIrsTHJRNWA56DeTbW1BxHR91M8PC7Hz+82Tp/+DA/3K86/ieUC5McffyQ/v3XHkJ2djclkajZmMBhYv349CQnntk5Q8PHx4ZprrumwLV5eXmzZsgWtVkt6ejq333478fHxQHMp42+//RY/Pz82bdoEcOZDqClKmqYMvb4AWTai0bhhY+PdOzIYnUSttmXUqCUcOfIIySkvYJYNBPi3bODR8kINXP2W0r3q+8cbulctO7sbu7vZv1TZxR3dUtKgvr6YgsIf8fe/rcsfpkd/3ojKbEaeOAWTTofLkCASs/VMzjZw6ZXtl8HuTURE348ZPuxFbG0DSUp+BqPxAmgA0Y2c6+TbGu8KBoOBBx54gMjISG6++WaSkpLOHGsqZRwZGcnWrVt5/vnn2bFjB87OzRVKTaZaamtPoNPloFJZY28/DFtb/37t5BtRq7WMGvUxHu5XkJr6N7Kyvmj/xZE3wf3bwNoBvpgFez4626azu9BXKbrzEZabi+TmrkGW67usa2MyGkndswvPylrSrK0ZNGgQS46eRKeWiDmtY0h01yp5eor+/xd2EWNlZU/4yH8Sn3AraelvMDLsnb42qddoK/JevHixxYjZ2dmZefO6Vxl78eLFeHt7c/jwYcxmM1rt2XZwTaWMQ0JCSEhI4IcffuCvf/0r06ZNY8GCBZjNRvT6fAyGMiTJCq02AI3GpU/1yTuDSmVDZOQHHD32JGnpC5FlI4MGtdSRsYj3SKXefu0j8NMLkJOgSFRbEBvrFImrlF7DFpuLmMjJ/QpX10uwt++amF7m4QPoDfW4BA7leGUlky6/nHfzjESUGrlsbCAqdf+MnfunVYIzODvHEDT4YfLy1lBY9HNfm9NviIuLU3YiNkGj0RAXF9ftc1VUVODr64tKpeLLL79s9VtDbm4udnZ2zJ07l2eeeYYDBw5QX1+s7Go1lGNt7YGDQwjW1q4Dzsk3olJZExH+H7y8ZpB+/O9kZn7c/ou1znDrCrjib3BkDXx2pdIas6vIspK28RmlpInOobjkV3S6nC6XVAIc3bAWjdFE2chw7OzsWJtbSIGthonHdYRP7l8llU0Rjn4AEBz8OI6O4aSkvIxePwDr63uAUaNGMWvWrDPpEWdnZ2bNmtUjVTfz58/niy++YMKECaSlpbXakOTIkSOMGzeO6Oho3njjdf7yl7vR6fJQq22xtx+GVuuLJA3MJuNNUak0hI9cjLf3bE5k/B8ZJ9/ryMWKhMbcb5XF2SVTIeWHrhl0eo/SXGTs/a02F7Gx9sbDo2tBgL62lpMpR/HQm8moqiI6Opo/9HZ46sxc4uWM1qEXqoo6iSR3d66sE4wZM0ZuXNwSWKa6Jp39+2fj5nopo0YtGbAR4flITk4mLCysr83oEmZzfUOaRhEfs7HxxcrKqdv/vfrD70qWTSQnv0Be/ncEBT3KkOCnOvY+y0/DqruUNpiTn1G0c1Sd+CD89n5I2wxPJ7dIBdXWZrJ7TxzBwX9mSPDjHb93E45s+ZHNn31AQHA4yVpbhlx1Jc/VO3Bzch2vzonE3d+hS/fvDJIkJciy3GZjABHRDxAc7IczdOhzFJf8Qm5u62V7gr6hUXysuiYdg7ESmwbxMY3G+YL8UAalmU5Y2Dv4+d5CZuYHnDjxf3QocHQZBPf9rHQa27EIVt4EtaUdM6K6CI6tg6jbLOb7c3K+QpKs8O+irg3A0Y3rsNUbyHZxZciQIazNrkFjlomppU+cfEcQjn4AERhwD66uE0k//ia1taf62hxBAwZD5RnxMSu1Aw72Idj0oPhYf0KSVISGvom//52cOv0J6cf/3jFnr9HCnPeVhdnMnfDJZUrns/Zy8EtFMXVsy0Vhk6mO3Lw1eHpOw8amawJj1aUl5Obn4OTgSpVOh9+IERx282RSroHLruzaAm9vMGD/Eis2bCD9ijiSw0aSfkUcFRs29LVJPY4kqRgZ9g6SpCYp+RlkuftLCQXt54z4WN0pQGoiPta/tr/3NJKkYkTIawQG3EtW1n9JS3sNuTUl09aIvRfu+0lRQF06HQ582fY1ZhMkfN7QXGREi8MFBZswGiu6ZRH22I+Kf6kLGYmDgwPfnC6m1kpFzCkdQaP6Z0llUwako6/YsIGC99dhE/E4DrM/xibicQreX3dROHut1o8RIQupqDjAqVNL+tqcixJlV2s+NTVNxceG9br4WH9CkiSGD3+ZQYPuJzvnS1JS/9ZxZ+8fCw/9pkgofP8YbHhSESlrjeNblTx/a7o2OV9ibz8cF5eubzY8tm0z9gaZHLOZqOhoDtp6Elpu4qpxg1Cp+n9qbkDW0ect3YzdqDtQqZXISbJzRzvqDvKWrsN51qw+tq7n8faeRVHxFjJO/ht39yk4Oob3tUkXBUqv1gp0+nxkswGNxgUbGx9Uqv5bbQHw+7oUHPcW4i5DiQRV472Ycl1ot88jSRLDhr6AStKQeeojZNlIWOjfO1ZpZO8Bc79T5I53Loa8RKV7lSWp6/2tNxeprEqkquooI0Je6/IaSeGpk5TVVOISEIwkSeyoM5DjYM19CdWEPeLfpXv3FgMyorcOuvKMk29EpbbGJmgaxR9/TP3p031kWe+gCJ+9jkbjxrGkpzGZzhP1CDpNvaGc6uoUKiuPUFWVTE1NOnV1WUiSGju7IdjaBg4IJ++3pxBPWUKFhKcs4benkN/X9YwyqiRJDBnyNMHBT5KXt4akpOcwm40du4naCq58Vam5L05XpBMyfmt+TlkmpG+GmLstNpvPzv4Stdoen25ox3n021UgQ7GnN8OHD2eP3h5XvZkpni5o7fv3v38jAzKi19i4Wh7XulBzsJyKH/6NUVuF+9UTcZ45A423dy9b2PNoNC6MDHuHQ4fncSJjESHDX+prky4o6g3l6HQ5Z7pmybKxQZvGFa3Wf8BU0jjuLWS7j4YPQmwo0Ep462QeTdMzZm8h9EBUDw3OPvgJVJIVJzL+iSwbGTlyUcc/FMNmgWcorJoLX14Hca8ofYm3LYSKLOUcu5b58fr6UgoLN+Hre0uX02my2Uxqwl7sbOwpMJmo9/bmiOTIdWl6xs1uuS7QXxmQjr7WWIm9xrnFeL2pjlInZ3xcr8cGqD1QQMWm96i3q8Hr6nE4z7gaK1fLHxIDEXf3KQT430VW1n/xcL8cN7ceEovqp+TlryfjxCJ0+jy0Nr4MGfpMtzVUr9fnI5tNyLKMqokSoslUPWCcPMB+byv+HqFFp1ZszreVeDNCy4vIjO7huYOC5iNJVhw/8Q5m2UhE+OKOL1R7DFd0ctY/CltfUZrNNy1C2Paqom3TRBo7L28NZnN9+4TX2uD0oQRqjQY0Q8NwdnZmU5EBlSeMqZFx8+tfSqPnY0CmbvZV/8JGb4lrp9gzdpoD106xZ6O3xL7SzfyavYT1pz8gvmYnObZqNMOuwtH/Bmr2O3LqkQ85dt9LFP9vLabq6r5+G93CsGHPY2cXTFLycxgMlX1tTq+Rl7+elJSX0OlzARmdPpeUlJfIy1/fpftmZmYSFhbGn//8CpMn38pjj73CZZfdxvjx1/P3v3+AuaH5eVBQEK+88goxMTFERkaSkqKkQn777Teio6OJjo5m9OjRVFUpYnT/+Mc/iIyMJCoqihdeeAGATz/9lLFjxxIVFcWNN95IbW0tAAUFBVx//fVERUURFRV1Rha5I5jNZvZuOcGHITZnnHwjOrXE+yG9Uxk0ePCDDB/+MkVFP3Hk6OOYzZ1IM9o4KKqXWpfmTh7AUKdE+A3IsonsnK9wcRmHg0PXI+7ENd+A2ppSa2uCI0Zy0M2LS/IMxE3rf5rz52NARvRfjLcmNVCLoaEnY2OUEuxgwsnBiPWJUwTnFuF7eifWKi0+3lH4Wo8kIPASnNUa6vbUUfn9p9TY1uJ5TTRe18ShaiJUNZBQq20ZOfKfJCTcTFraa4SH/7OvTeoW0tJep6o6udXjFRUHkeX6ZmNmcx3JyS+0uqHM0SGMkJC/tTl3amoq77//Ev/618uUllbg5uaMyWRi9uwHSErKYMKESAA8PDw4cOAAH374IYsWLeKzzz5j0aJFfPDBB0yaNInq6mq0Wi0//vgj69atY+/evdjZ2VFaqmwKuuGGG3jggQcAePnll1m6dCmPP/44TzzxBJdddhlr167FZDJR3cGgJHFfDqU/ZDBMB/nTLG/kKdT2nhTDoMB5qCQNqWmvkHhkPpERH6JWd7CZjiSBrqWIHdCs+XxJye/odFkMG/psFyxWMNbXk3EiDY27HyqVih9LDVR7qYjN0jP4Dvcu3783abejl5Sl83ggR5blayVJCga+AdyAA8BdsizXS5JkAywHYoES4FZZljO70+i8YTdhMDU33WBlRfGIO/ltciynq06zK2cX+9J3UHIwicDsHWTl7WO/So2HXzgBmnAGeYfjamWL4Q8DmRuWUaWtxX16JP7XXIZkPbDqoJ2doggKeoyTJ9/FwzMOb68ZfW1Sj3Ouk29rvL2YTHoGDfJj7LjRgMzatT+zbNkaTCYT+fnFHKTMc8IAACAASURBVD9ezIQJyrk33KB0W4qNjeW7774DYNKkSfzlL3/hzjvv5IYbbiAgIICtW7cyb9487OzsAHBzU2R0jx49yssvv0x5eTnV1dVMnz4dgF9++YXly5cDoFarW8gdt0ZGajEnvk0lrNLMKUcdDwRXg2Q5R+3c8O2htwgImIuk0pCS8hKJiQ8yatTHqNW2bV/YFOeAs7n5c8cbyM5ZgbW1J56eV3XRYkjb+hMGCfRevoSMGMEmWy+GVpq4emwg0gAoqWxKRyL6J4FkwKnh9TvAYlmWv5Ek6WPgT8BHDc9lsiwPkyTptobzur7/uAmlJstml5qsmJ6QxmhHO0a7TOeRuBsImqEisfgwOxM3k7V3PxUZCRTVJpKgAvvAYQxRjyTYJRB3ayfkP8yc/OFrKqyrcYqLYMisS5E60R6uLwga/AglxdtJSfkbLs6x2NgM7AXotiLvXbsmN6RtmqO18SM25qtOzWk216PTZWNnZ4u93VCOn0jl/feX88svX+Hu7s78+QsxGs/+PdjYKFGpWq3GaFQqS1544QVmzpzJDz/8wIQJE9i6dSuyLFvM6997772sW7eOqKgoli1bxq+//topuwtyKzn0TTKhhXpkexNPhFTxR/Ag7OvsCS7I4rSHLyZ1k/8zssyEU8fJz4nEx9+p9Rt3M/5+t6KSrEhKfp7Dh+8nKupT1Gq79t8gbgFseEJJ1zTSpPl8Xd1pSkp+IzjosW7ZtHZ043qwc8IAJGvsOOVozd0Hqwl/pP+qVLZGu3L0kiQFADOBzxpeS8AVwJqGU74Armv4eU7DaxqOx0ndvHrlb2N59d5RrcJJrebbgjL+nJLFZftSCP8jhbfzXKkZOo+4h5dx87++IPDh69GHu1Oen8axk9/zbeH7LJf/y0HpADpbezztR2GzR0XmE9+S8ORnJK3Z2SMNLboTlUrDyJH/5P/bO+/wKKvsj3/u9Jn0ShLSKKEkpBI6qHQBFdfeQd3FsvbGrlv0t65l1V3brnVxYSmKDUXERaUovaSQAKGFlgQSUkidPnN/f2RAkAAJScgkvJ/nmSfvvPXcvDPfue+5557jdlvZUTCzZdPQOyE9ez2BSnVqj1ClMtKz1xPndT6324HZvB8p3ahUOtRqA1YL+PoGEhMzjIaGIJYtW37O8xQWFpKcnMzMmTPJzMxk586dTJgwgQ8//PCED/6466auro7IyEgcDgfz588/cY6xY8fyzjvvAI2FVGprmx57qa22sOzdLOrfyEbbYOOphHruGBlEdlQ4w3fmctO6b5i4M4vLduXgazWDlBjsVhCCw34+bJ+VjdvdwklNrSQy8lqSEv/OsepN5ObehdPZArfUOYrPl5R8hBAqotogr425tobiilJkWBTBwcFsdvrjb3czOiwQnbHzebyba/HrwFPA8efAEKBaSnk8QLYYOD5zoDtQBCCldAohajz7V7SJxcCdbiMvOe04ND//fmidkke0vvw2vSduKdlrtpFTayanzkxObQPvFZXj8IhfN90Q0qeOJsVHj75oG+51S2HbPvY4l5NrcFIXqSfTnUaSqjfhhgTEFknRmqVUiGrUw3qRNHUwOr333Wwfn5707v07du9+lpKSBURHt66ajjdzPLqmLaJufhZ5JwZDNND4uUpNTSU9PZ2kpCR69uzJiBEjznmu119/nZUrV6JWq0lMTGTSpEno9Xpyc3PJzMxEp9MxefJkXnjhBZ577jmGDBlCXFwcycnJJwZu33jjDWbMmMGsWbNQq9W88847DBs27CR7JdZ6O6UvrMXgZ+RPfe2s6BGOzmEgc89WksuK0Luc6NDjRJJQXkJCecmJ41f3TiE/pjfHysys/Gg7Y29NbvH/rDVERExFCA3bdzxKbu500tL+0/wwyDMUn3e5bBw+8imhoeMwGCJbbeO2hQtwGoyYjT74x/Vgq68/kwptDJvaPiGp7c050xQLIa4AJksp7xdCXAY8AdwJrJdS9vbsEwMslVImCyG2AxOllMWebYXAYCll5S/OOwOYARAbGzvw4MHmJ+ma8/Rafgqq5sekYOoNBnytVi7dXsWIMn+mPpaOf6gRre5Ul4vV5WZHvYXsOjO5tWZy68zsNf8cAdDLoCPeVotxTy4B+bmElx+hzmShMkwyxNmPjIYkgnwTEGodToeZclmNOzOe5KvSMfl6j09fSknu1juprt7MkMFLMJl6dLRJzaYjUu+63U7M5n24pQOTMd7rinOfjJSShsp61BYXO4sP8PKuKlbFBqFxOUkuLiSlpBCTxUGIIZrLJo4mcVBP8vPz+XLRItwnfc8dKjWL0i4BjZH/rrcQfFsSPfte+HwtR8uXsW3bw/j59ictbTbaJkKmm8uRI1+wo+BJ0tPmtkmY8ZxpN3DUNwhrRDTbe2fyU7cIns2zMeORoa0+d1vS3DTFzRH6F4HbASdgoNFHvwiYCER4eu3DgGellBOFEMs8y+uFEBqgFAiTZ7lQS/PR//2Rj6jz3wOqkx473Sr8ahMwWBt900Z/Hf4hBvxDjSf++oUa8A8x4husR61WUeNwsrXO4un5N5Bda+aovfEhRSMlkbUVhBYVEnG0GIPtABZtEUMaYhnWkEawXxIqnQ9ul4Oj7mosKd1JmjyA4LCOFwqbrYwNGydhMsUzMOOTTlGTFC680LvdTsyW/bjdNo/Ie2eqWel2YymvQWWTODRayvSSfYf2c3eFnQEl+0g7uJeAOiO9e2RwyRUZhMed6nfPy8tj+fLl1NTUEBAQgLnqGEf8A1mUfikjyh08UGBm6B8vQae78J+T8orl5Oc/gK9PAunpc9Bqz2+ey+Yt1+B01jN0yLJWz3OoOrCfD3/3MJY+acQlJvF8UAL9qly8ldKT2CTvirZprtCf885KKX8P/N5z0suAJ6SUtwohPgWuozHyZhpwPIB5sef9es/2FWcT+fPBHHAQxC98iyo3dQG78endgE6YEC4dtRYdVftV2HLUCKcW4RmSEAJ8gwz4hxrwCzUyPMTA5aEm/MKDMAdo2Skd5NZbyKn1IyconBx3o/k6m5Ud5SV8ZTmMtHzFpcVmJpX3Itw/FVWBnYYdW9jvqqWmXzf6jO9DdFxgWza72ej13ejX9zm2bX+IAwffaXXBha6IlC4slgO4XTZMpjivFHm3zYatohrhUOHQ6TnqA2atGpDonQ6mr1xDiCWa1EHXkDE+noCwpqNYUlJSTqm8VVFayttvvcmgfdv5KSGFSyrA/GEul997Tr1oc8JCx5Ka8i55+feSnXMb6Wn/RadrmZjW1uZRW7uVPgl/bpPJbHnzZuPwD8alUrHObaBGpyKzqIGYm04vOt5ZaM1P+EzgYyHEX4EcYJZn/SxgrhBiL1AF3NQ6E0/HJaxNbxASlc5NZU0JDQ0Njes0QGjjtGyT0QejzgeNMGJzGThq0XK4QI2jTo3KpUclG/8dao2KHiEGUkMNPBBioCZUxyE/FdvVNeSFRZCnjcetGsXWFPjQWodfw376FhcwsVjHYHs03fbaYG8+Oa46KnuEEHtJD3onhp4yw7K96dZtChUVyzlw4C1CQy7F37/tS+x1VqR0YTYfwOWyYjTGelXWSel2466rw1ZVjRQGLHotlX5qrFoNIDHZbfhY3djsaq4YOpXky6Ix+bfMdRgaEcHooSNxbt5AUXAEr/YNZu4GC/mbS0gedOGTdIWEXEpqyr/ZmjeD7JxbSU+bi14f1uzji0sWoFIZiYy8ptW2SCnZtS0XGZtAaFgY3/tEEFvn4qrM2E41I/qXtEjopZSrgFWe5X3Aafk/pZRW4Po2sO2MBAQEUFNz+uSJgIAAZsyYAYDD4aCmpubEq7q6+pT3NfWHGiMO9DS+AJ1Wj1Hng05los6lp6Zai7NIg9uiQe0ykOzWkoIAk5Oq8FoOmKo5YBQcDo1iVWoKq1IBKQmtLSPtSB3pdQYGHHOgnb+L7a5sjnQPpNvgaJIGRaHRtL/o9+nzLMeqN7J9x+MMHrS45XHLXRAp3ZjNB3G5zBiNsWi1Fy688Gy47XYcVVVY6+pw6X2xGUzU6LTYtDpA4mN34GMGrUqLyV9LdZCexBE9z/t6o66czI6tuYzek8NnGaP5c7KOF77aSUP/UHx8WziZqQ0IDh5BWuostub9huycW8hIn9esEGGHo5qyssVERPyqTX6wSzaso0avx27ypTy4G4X+em7Ja6D/PZ0jS+WZ6BzO218wduxYvv76axwOx4l1Wq2WsWPHnvI+NDSU0NCmB5ncbjf19fVN/ghUV1dzrKYcm8sGJhpfgEql9vwQGAl06Umr0pFUJ5DmnZhdJZT52ykNi6S0Wzzr42P5wdD4hdG43PQ5piK5HpI2FVH23S4cJg2h6TGkjozhufeWsqR3NEeNKsItbq7YW8zzD7c+Z4tW609i/5fJyb2DvYV/o2+fZ1t9zs6MlG4sloO4XA0YjTGtGvxrG3skztpaLFWVWKUEnRGbrx8NOsNJAu/C1yzRatSYgnUYfLQIIdqkd3nn44/y6p//zMg9ufyQNJgl8ZJj729h8mPnji5qD4KChpKW+h9yt95FVvbNZKTPw2CIOusxR458jttta5PiIgBbP1mAIzAUjVrNBlUIPg7JuLAAdIZOKZUn6LTFwX85wDR27NhT/JBtgdVqbfJH4Pjy8XC4k9EKDdituC21NGh0lAdEURYcR1lwBKVBRhzaxp68n8NNYo0bndPN+nANzpPcOjqXm9u3F7WJ2ENjOoGi4tmkpc4mJGRUm5yzPWjPwdhGkT+E01mHwdAdna7j/K0uux1LRTkWqwWnRotUa3Cp1Jh1BqzaxjkiPjaJj8WNQafG5K9DZ9ScIu5t9b8q3LqVeQs/YvmAYeztFssHm82EDopkyNjzf1poLTU1OeTkTkerDSIjfR5GY9MTlKR0s37DWHS6MDIHftLq6zrtdv51+3VU900jqG9/Xo7ox7j9dt6YmnLG8Y+Ops0GY72VXw4wtQcGg4GIiAgiIiKa3O50OqmtrW3yR+BYVRWypoYYRy0xZflQlo8bqDH5U+kXRplPIAf8gygL8WkcHT4Ju1rFlwlRPN9G7ejV60kqq9ZQUDCTIUOWotV2zCBxRyGlxGIp8oh81FlF3tfXt8W5ZZplg9uN9VgV5rpa7EKFVGvAYMIlVFh1RsxaLSBPCLzRoMEUbmj3nmSv1FTS1m/AtmcrpQGh/HGAlrdXHaAypRshHRRBFhCQTkb6XHJyp5GdfTMZGfMxGmNP26+qag0WyyF69ni0Ta67Z9FnmINCkCoV66U/TpVgsNnttSLfEjqt0HsDGo2G4ODgE7lLfomUkoaGBooK97JzyyYO7t6JvqqMMF0Ryf4BuFQa3hw+pcljK9vwC65WG0hK+jtbtlzLrl3PMGDAG2127o7k89IqXtx3hBKbg+56Lb/vGcm1EafeCyklVmsRTmctekNkiyM6WoOUEntDA+ZjVdhcLtxqDeg8yfOEGovORL0nisZkc+NnlZiMGkwRRjS6C5d648rf/Jr9M2cyelcWi1NH8u++Ptzw/hbG/37UBQ0gOBl//xQy0ueRkzvthBvnl3NCikvmo9WGEB4+sU2uue27pTjDoggOCycvJIKMow6umtB5cs6fjU6ZprizIITA19eX/qlp/OruGTz80t+5++HHGJmSSEBJIfr8DfjaLE0eq3G7yN5a0uS288HfbwA9ejxE2dEllJYubrPzdhSfl1bxxK4iim0OJFBsc/DEriI+L606sU+jyJfgcNSg10eg1zV/UlB9fT1jx449kYb4q68ao4cbGhqYMmUKqampDBgwgIULGzNl/u53vyMxMZGUlBQee+QRakqPkLV2LWMnTmTEhIlcd9vtHC4tRSV02PX+lPv4Ua9VYbK56VbrJkKtISzSB//QCyvyACq1mtsfe4SY4n2kH9rNN931HDHp+GnRrgtqxy/x80siPX0ebredrOxbaGjYe2KbxVJCRcUKukfdgErV+sFja0UFBxwOXAYTBX7hVBrUZBZZ6d63a9Sv6LQ++s6O2+2iaFs+f/pxEyv6Z+A8KemUyu3GLQQRZif31Eruu6X1xY0br+kkO/smGsyFDBm8tE2mirclJ/ud/7SnmG31Tf8IAmTVmLE38dnVCcHAABPIxiRlUjpQqXQIlZYBvkaeSzh7Qqrjrhun04nZbMbf35+KigqGDh3Knj17+OKLL/jf//7HBx98AEBNTQ0ul4thQ4eyed1arA0NVNXW4h8cwrTp07liyhRuv/U25sz7iG+/X8YrHzf6kk02N/52ia+vDqOvFpW6ZX2u9hjP2LB4MUs3rOfLweMxG/2Zvb6euGmpxPToWLGrb9hDTs5tSCmJjbmLkpL5JxLaJST8idiY6a2+xubXXmHZoUOIoDC+S55AnVbHbN9QUi473WXkTTTXR6/06DsIlUpNXEoa/fPWcOnOn5NO+VrNjC7IYmrWKuo0Tv7aTcNj76+ktuYMcwdadE0NiYmvIqWTHQVPIeWFTWjVljQl8ifWS3DLRpEXKi1CtLyup5SSp59+mpSUFMaNG0dJSQllZWUkJyfzww8/MHPmTFatWolGurHVHEOjUXPPAw+y+PsfMPj6oVeryc7O5rrrb6VOGBh1xx1s3rgBk81NRL2baIOO8EhffAL0LRb59mLIlVfSQ61mTMFmbCp4cYCR3XOycTo79nPi65NARvpHuN0OCve9ckrW0sLCV1tdbAZge9ZGnH5BWGJ6sivQwIj9FvoP79whlSej+Og7mBCdmgFbV9Pn8D6kVodw2NGVl2C0NZBqszGnbxoLEiLY820Wv+vdnRGZ8a26nskUT0Lvp9m5648UF/+XmDboDbUH5+p5Z67bTrHNcdr6aL2WjxJ9sdsr0Om6oddHnFco4vz58yk9cpj/ffEZKiEYfOloqivK6dM/kdUrV/L14sU89eRTXHrpJTz62OMsWfING9as4eulS5k/fz4LP12EBI4YdGi1KnzqHaiBGJMevUnjlZNvhBDc8PjjlP/hdwwPCufHvulsipA45+Yy8c6MDrXNx6cnarUBl+vUTJ5ut4V9ha+2qoRk1fbtHPH1B5WKPEMkBqdkYljAafmyOjPe0ZW4iBl10x2YrA34FubjtzML38J8DA21uOx2NId38ZL5CKN2byMrzMhvjh7ltbmrW51aNirqJkJCRrO38OVT/J6did/3jMT4i+IPRpXgyVgTdnsFWl3weYs8QEVZGf5+fmAwsTorh6KSEmqrq9m6ZTN2h4Op117Lvffdx7b8bWCzodfpuOm223n2+ZfIyc2lVKcjZchQVn38Cd2tsHbJF4y6ZNSJOHhvxScwiKvvmM6AXdnElx/mnwk6ONzArryyjjYNu728yfVW25FWnTdv7n+wB4WjDwgmNyyYkUV2Lp3Qu1Xn9DYUoe9g+o8azYQZD+AXGgZC4BcaxuX3PcxtL7yGWqMh53+LeMDHwbQNa3Dg4OUoHx6ctYLyo+cfAiiEoH+/F1GrTWzf8Rhud+uqMnUE10YE82rfGKL1WgSNPfmXegUwyb8erTYIgz6qVYI6eeJ4tuZvY9KUKSz68kt69+6N1GgpKNzHFVddxcQJE3j7n//kL3/9K+j1XDF1Kv2Tk5ly1RQef/FvmOxuXvvrq3zz6QLGjBvKRwsX8MYbnSPaqc+QYaTE9WD0tvXoHTb+lKzn8Of5WK2nP0FdSAz6pseUzrS+OUi3m22Hi5F6A/khcdjVgiFmN37BnbO06JlQBmO9GJvZzPcf/JNd634iNjkNu7E773SPY2+3KJIq6ngiNIxJo88/P/bRo8vI33Y/8fG/pVfPx9rQ8vOjNQOMNls5NlspWm0gBkN0q0ReSsmRkhJoKrRQuonoFoFK0+j1rKyro9YuqdNqkICP3U2QUBPgr0Ojbb9H//bO9Gm3Wnj38UfY2juZb9JGcv1BO5PqXEx+cNi5D24njheEd7t/HqRXqYz06/f8ebtuipYuYc4P3yH9Q/h08BVENEhmp/UmKqFzzDVRBmO7AHqTiSkPPcn4GQ9yeFcBtbvW8orawvjtO9gZZOJBWw3P//t7HPbzq34VHj6RyIhrOHDgHWpqctrY+guH3V6JzVaKRhvQapF32u1UHC45bRLbzwhUGg3VDQ3sq6ylRKqp1WrwsbuJcQl6BJkICTG2q8hfCHQGI9c+/AgJO7NILi7k0zgdtTZJ9ppDHWZTZMRU+vV7HoM+ChAY9FGtEnmA7K+/xOkXyJGoXhw1qhlUbCGyd8emxmgPFKH3coQQpIydyK0v/AOjnz/rPp/N3X5m7lmXg8Zp562eocyY+z0HDhw7r/P36fNnDIZItu94HJfrwhaMbgvs9iqs1sNoNP4YWyHyUkpqqyo5erQMhzjz10Ii2FdZyyGXoM4j8LFuQY9gH4KDjKi9JIKmLejepz/DLhnNiJyfCK6v4ZkkLdXf76G2+sxhr+1NZMRURoxYzdgxexkxYnWrRN5VV8du6QahoiCkB+EWNzdmxnv1GMr50nU+lV2c0Jg4bn3hHySPHk/W0i9JcOzihX3FJB05wrc9I7g5O49Pvmq5+0uj8SOx/ytYLIfYs/fFdrC8ZbTElWi3H8NqLUGj8cVojEGcRaDPhs1ioaykhHqrDYTAZLfjp1Jh1Wip9PGn3DeQSh9/agw+VPoFNAq8o1Hge4b4EBRoRKW6cOJwId2tI2+8hVghGZe3jjqdin8m+rL6/U0X7Prtyc75czGHdMPqE0RBsA8j91tIHHr2JGqdFUXoOxFavYEJ9zzElIeepKLoAEdyv+KvGidX5BZS5O/DTKObP769GIu5ZYNmQUFDiI29m5KSBVRUrGwn68+NwWCgsrKyWULmcFRjtRaj1vhgNMadl8i7XC4qSkuprKrCLQQ6u41Qg5GAuDiqtT7UG0y4Ped1CxV2jRa9E2Lcgp7BjQJ/oXt/UkoqKysxGC7MYKFao+Hqx2cSuzePIYXbWReuZ49Jx5pvdl+Q67cnm7dsQmr17Izqh87VGFLZ2V1uZ0IZjO2kHCs9zDdvvEzZvr0MGDOZrCOhfJQRQbmfP5ccLOGPA9NISW7+hA+328bmzb/C7qhiyOClHZLd0eFwUFxcjNV69slhLpcFh6MKlUqHVhvSYpGXUmK1WLDZG6ONVG4XRrUGjZ8fwjMAe9hix9WEhqslRBk7tkawwWAgOjoarbblE8HOl5xlS/nf5wv57PJbOeofwgfra0n/dSbdIr2naEtLMO/fz2tvvY41IJS5w69gyGE7s65Kwzfowufibw1tVjP2QqAI/fnhdDhYvWA22Uu/Ijy+F4awMfwn1EB2XDcia2u4t8bBb24b0+zEVHV1BWze8itCQ8eQPOBfXumrrKhYSV7+ffj5DSA9bXaLSgBKKcnfupWlixZhFQJtXTWD7XYunfk0uujGCVrL9mzjvc1lrIsIbnJAVkjJkTHpbdaezoKUks+ef4asehsfjbmWaLObP26rZ+LTzf98eRM/PfNHVqBmX3QS3/Xuw+8KLDxyf8dFFJ0vStTNRYBGq2X0tN9w9VN/orbiKEfy5/Goo57rNpVQaTTxfGQAT7z9GdVVDc06n59ff3r1fJTy8mWUli5qZ+tbTlXVWvK33Y+vbx/SUj9skcgXFRXx9iuv8MWXX2K3WYg/uJf7rrmO8W+9jS46mi+2Z3HVxyu465CDTd2CMZ4hkCnc2nnTRrQGIQSTH3iUqLJDXLZ9M4UBBn6INLFiYX5Hm9ZipNtNblkZEsiP6kn/Kic3XZ7Y0Wa1K4rQdwF6DRzC7X97k/D4HuT98B/GqQq4b+VRQustLEjqw63frODHVduada7Y2F8TGDCIXbv/D4ul7bJntpZjxzaxNW8GJmMP0tPmNLsEYGVlJQs+nMWsWbOoqKnGv3gf1/SIZ9qs/xI8YiRzstcyZeFKHihTkRsaxIRDVl6rUXFNSRUG56lPuwanZNzR6nZoXefAJzCIyfc/TFL2SvodPsj8HkaqDtWwf09lR5vWIspX/cixkDAqAiI4YtIwuNhCtx7eUVKyvVBcN10It8vF+s8WsGHRJwRFdEcTeDlfdvdjXe8QQuprmH7oKI/OuOacA04WSxEbN03Bz28AGenzzjuapa2oqckmJ3c6en0EAzMWoGtGumGz2cyq779nc3Y2uN1oK0tJNuiZ8Ofn0AUE8s66n1h8WLI1JACDC8YdsnJFkB9XXZ2E8ETQPD57OT+EB3LUoCLc6mbc0Wr+Pn3sOa7c9fn27dfZnL+NeVdOR42Wv2+pZfzvxnSa3DCfPfoQ2wKCWZM0hiO+vnzkH0bKiJiONuu8UHz0FzEH83JZ+s9XsZnNxAyYyhpHd74a7I9drWLyrt38aco4YuLDz3qOw4c/o2DnTBJ6P01s7N0XyPLTqa3NIzvndnS6EAZmfHTOgtEOh4ONGzbw04oV2N1utNUVBFdXMOXeh4kZOoy/r1jOsmM68kP88XFIJuy3Mjk0gMlT+6G+AAXbuwJ2i5nZjz9AfrdefDZqEqNKzUyv1zJ5xpCONu2cuBsaePlPf6A8JJL5wydy1U4z7/x6KGpt57z3io/+IiYuJY07Xn6L6P5JHMj5lMGOtfxmWRUx1RYWJyZyx7pNfPrJ8rOeIzLyWkJDx7G38FXq6zumAEVdXYGndmggGenzzirybrebvLw83vr73/lh+XJkTRU++7YzvFdPpr09i9lVtYxdtJZXVaEc8vPl+gILbx3T88a0QVx5baIi8i1AZzRxxSNPErdzI5mFO/kp0pdCi5VtWw6f++AOZsfHH2P1D2R/9yQ0bsmkUL9OK/ItQenRd2Gk282mxZ+zduFcfAJDEcZJrOwdxsp+AQRa6rlx1x6emnEzPr5N18S02yvYsHEyen04gzI/b5NKPs2lvmEP2dm3oFLpGJjxMUbjmR+t9+/fz7Kvv6a0qgp9fR2qysN0MxoY89ATvFdwgDUEsSvQRLDNzfg9ViYG+jH66r4Y/To2TLKzs+bjuaxd+hUfRCgwPQAAGRxJREFUXXMfx3wCeHNjFeMfuQyTr/f+X9976AEOhnZj4dAppJU5+O9VafgEdK6QypNRevQKCJWKIVdfz43P/g2hkjSUzufyA9u5fWUNbpWOf6emcs+8z9m8YXuTx+t0ofTv/yL19QXs2//mBbPbbN5PTs7tCKEmI33eGUW+vLyc+XPmMGfOHI6VlOBTXIi+eA+Dx09i+yXXcMuuemYFdqdKb+DmrWb+VqzmuZvTmXx7siLybcCw624mMqI7k1YtwY3klUQTP3ywvqPNOiMN+/ZR6h9ASUgcDVoVw8yuTi3yLUHp0V8kWOrr+O7dN9i7eQNBkUmUuy7j2+E+7Izwp+fREm6rquGeGTeiVp8+oFZQ8HsOH/mUgRkfExh4zs5D6+y0HCIr+2bcbjsZGQvw9Uk4bZ/6+npWLl9Odk4OaocD/yPF2M1VhMcmkJs4inXB3TngpyOywcWYXVZG6kyM+lUCodHND8dUaB6VJUXMfephsgZO4Pv0ody6p4obescxZEyvjjbtNL597i9scLlZMnASRqeaz9ISiOjRuROYNbdHr1SYukgw+vpx1eN/IHfZEn6cOwtfw2Gu/3EiWxK7sywpktf8Ail87X0evfFqomNOze+dkPAHjh3bwPYdTzBk8JIWxa+3BKv1MNk5t+NyWcnImH+ayNvtdtavX8+aH3/E6XTSrfgQVmsNQmOi9NJpfBwTT7GPhuh6J3dsamCg1DHq2kSiErpGgWdvJKR7DJfecReO/7zHvsgYPuoVSfqmvSQkRxAc5tPR5p1Aut1sq66mLKoHJb56btleR8SvOrfItwSlR38RUra/kCWvv0RNWSmmoFEcCh7EV8OMHDPqGLpvB9PCwrlq6vhTZsZWV28hK/tmoiKvo3//tk9+ZrOVkZV9M3Z7JRnp8/D3Tz6xze12s3XrVpYvW0a91UpESQmahirMKiPFg6byfUJPjhjVxNc6uHS7jWS7huFTe9EjLdQrZ/d2NaSUfP7iMxQUHmD2dfdhcLl5Lr+Bq2eO85pZs/uXfcec9etY2+8SigIDWRgYTurQs5er7Ay0mY9eCGEQQmwSQmwVQmwXQvyfZ/1sIcR+IUSu55XmWS+EEG8KIfYKIfKEEB1bbFLhNLr16MXtL71BvxGX0lD1E9Hli5j+bTkDjlSzrvcAnne6efHVt6mpqTlxTGBgJnFx93D4yCeUl//QpvbY7RVk59yO3V5BetrsU0R+7969vPuvf/HVV1+hKy4hcdtWXA7J7sG38sGtj/HflAR8bC6mr6njnvV2pl3Sk1v+NJie6WGKyF8ghBBcft8jBEkHl69bQaXJxIJowY9f7uho007w44rl1OqMbA8PYtghKwMGdc0slWeiOa4bGzBGSlkvhNACa4QQ33q2PSml/OwX+08CEjyvIcA7nr8KXoTOaGLSA48TOyCV5f95F711PlNWTqBHquSbxHBmJQdQNGsed40czqDBjbldevZ4iMrKHynY+XsCAtKaNXHpXDgcx8jJuQOr9TBpaf8hIKDxWqWlpXy/bBmF+/fj09BA+s4CjnVPYPOld7MkIYpjehV9q6xcucFOjzoYODGOlNHRaDrJpJ2uhm9QMBPveRDbP14gI7YPa3onMDj3IAkHo4mO69hqTZaqYxzSajkY0QchYUqob5eqG9Aczin0stG3c7xAqdbzOpu/ZyrwX89xG4QQgUKISCll6yr4KrQ5QggGjB5PZEI/lrz+EhVFXzAgdzBRlaNYlKHly7ThFG3fydSsrdw2/UaMRiNJia+yecvVFOz8AynJ77aq1+xw1JKTOw2zZT+pKf8mKHAQtbW1rFixgtzcXHQOB6n523B378+Pkx5kSe9warWCxEoz1653El3lJmV0NBkT4zD4XLhMjgpNkzBkOEmXjsOxcgGHwh/l7cRQouet5saZUzp0nsKmTz7BptOzNTqWzFIHU65K6jBbOopm/feFEGohRC5wFPheSrnRs+l5j3vmNSHE8Til7kDRSYcXe9YpeCkh0THc8sI/SBl3OfaGTQQd+Jjbv6tkYHEZW+L780637vz19X9RuHcfvr596dXzCSoqfuDIkV8+zDUfp7OO3K13Ul+/m+Tkd/DxyWTFihW8+frr5GVn03fnLtLK3Ky84mGeuuZGFvTvRlxNPfeuquW65VbG9g7jtr8MZfg1vRWR9yJGT59BSFAQVyxfhEOt4V99TXw7t2MLlWw5sI+DIdHUa9UMMzsvytDaZkXdSCldQJoQIhBYJIQYAPweKAV0wPvATOAvQFNdvNOeAIQQM4AZALGxsedlvELbodXpGf+bB4gdkMp3772Fqnou41eNp/cQDYt6BvNRxihKl/yPK2KimXLlbVRUrGD3nucIChqC0diy++dymcnd+mvq6raRlPgm+/f5sHL5a5itVmIPHiJShvP1+Hv5Js4fuwpSj1YxapeewCMOeqaFMeTengRHek9Eh8LP6E0mJj34OJ88+zsuyc1l5cCBfFtxkD75ZfRJPnv6ivbgYHY2tQYju6L7EV/r4q6JAy64Dd5Ai56npJTVwCrgcinlEdmIDfgPMNizWzFw8gyXaOC0udFSyvellJlSysywsLDzMl6h7ek7bBS3/+1NwmKjcTZ8TfzalTyYVU+w2cI3qcN5r97Ka//4J6GhTwCC7TuepLEf0DxcLitb82ZQU5NNUNCTfLpwF9988w0+xSWkFbvYeOldzLj9Rr6M9yf5cAmPrK3mylWCRF8T1z41kEn3Jisi7+VE90ti0NTrGLh5EfFHK/iydzQ/LFmN1dqyymdtwdqlSyn1C6bI38iQww2Ex3btLJVnojlRN2GenjxCCCMwDtgphIj0rBPA1cDxPLiLgTs80TdDgRrFP9+5COwWwc3PvcLAKVfjsm1FnT+HO3+qYdjhInJj+/BxvwzenPcd0n0zNTVbOHjo3806r9ttIz//Po4d20BlxZV89WUxtqJi+lRoyR58M7+9+Rq+iTExsHAPj6/ew6R1JuIwMeW3KVz9WDoRPS+euOfOzvDrb6Fbj15M/m4+OpeLtweEsej9FRfUBovZzF6bjX1R/fCzu7kro+cFvb430ZwefSSwUgiRB2ym0Ue/BJgvhMgH8oFQ4K+e/ZcC+4C9wAfA/W1utUK7o9ZoueyOX/Ormc+g1VuwHZnN6BWH+HVpFdUmXz7LHMvsbYHU1fZh375/UFdXcNbzud12snPuo7LqJ3bvHkLhVhNRDb7kpl3DE9dM5vsoA0N25vPIsh8Yl+VPmD2acdP7c8MfBhOfrMTDdzbUGi2TH3gCf1sNl69eRYVvAAsDrGxZs++C2ZC1ZAm1Rh+2h4cxtMhKcubFFVJ5Ms2JuskDTqudJqUcc4b9JfDb1pum4A30zBjEtFfeYskbr3B417eELD3IH0b+ivfDGlg2YChlh0O43/ACWVn3M3Lkt2g0pxetNpvrWLf+LoTIZt/uQVhKhpAzYAg/RegxOF0Mz89iYEEWvo5hmAKHkXlDPAMu6X5RZBXsyoREx3DJrXfimv0eqfGJbOgZy1eb1tNnQCT+gU0n0msrpJRszs+jMD4VKWByqC+qiyyk8mSUmbEKzcLtcrH+84/Y8PlChCqQoF7X8780+DEklihzJY8bn8GvvieRkQ+yZs0aampqCAgIIC4uBqfzfUJCC8ne/yvW+F3PxjA9RoeDwQU5JGf9gJ9IQu8/nPTxvUgbH4veqGTm6CpIt5vPX3yGwp27mHPDo9i08GTeIe596sZ2ve6BXbuY9fFCFg6ZTO8qF19MzeiS0VlK4RGFduHQtjy+fu1lrPV1GALHYLl2EO+6NTi1krt4m0OHklkVMZx6vRFfm4UrLV/SI3AXiy33ss0UiclhZ0RxAf1XLsbo9kPrM5Hk0RlkTo6/aDIJXmzUV1Uy58kHKPKN5b9X3UjPilL+oAlh0rXtN4/y49dfZ4nWjx8SB/LbfQ386e4R7XatjkQReoV2w1xTzdevvUJxwVZU2gTir5vGy7KM/YHdUEkXbvHz7FQh3Uihws9uZZytlF5ffY6or0FjGEKfEVMYdnUfAsNNHdgahQvB7g1r+Pq1l9g45BZ+Sk9kSsE2nrvyCqJi2n7WrNls5tWXXuLrjAm4VXr+l96XsJiuGW2j5KNXaDdMAYHc8OfnGHHjNNzOQvYvfJnH9rrQu6yniDyAFCoMTgtP/7SU3vP/jcqiIyb1fm7+y4NMmpGmiPxFQp+hI0m8ZAyDNy8ktqKO7xL6MnvB57jd7ja/VtbKlRwJCKYowJchhxu6rMi3BEXoFc4LoVIx9Jrruekvf0Nv0nBoy7vYzlCByqo2ULV3KwGR47n26Ze44enLCY9TvnwXG2PuvIeA4BCuWL4AtRQs7BvP3A+WtOk1pJRs2ZLFnsg+mBxufpPpfXnxOwJF6BVaRfc+/bn7jbfp3j+DEMqb3CeEcsb8+jnu/sdDxA1QJsddrOhNPkx64DF8qouZsHYjZQEhfK41szV7b5tdY/++fRzRaCnoFsmQQzWkZnT+VMRtgSL0Cq3G4OvLjc88w/XWz9BJ6ynbdNLK9dbPSB+fglApsfAXO9H9BzDoqmvpu/Nbkg9UsDm+D7NX/oC5wXrug5vBxu++Y1dEPC6VYHJIICrlMwcoQq/QRgghiFlTwN2u9wiVR0G6CZVHudv1HjFrzj6ZSuHiYvj1txIW14Nxq+cQYHHxTf8BvPPOvFaft66ujoKyMnZG9SK5rI4bLtK8Nk2hCL1Cm1F7MIbolXt5pe5B5snreaXuQaJX7qX2YNPFvRUuTjRaLZMffAK9vZ4rV35LncGHr6PCWLp4VavOm71xI/vCoqgx6BhyrA69qevFzZ8vitArtBmp42+gujCEHQsS2Pp+f3YsSKC6MITU8Td0tGkKXkZoTByjbp5ORNEmRmw7yM7IOD45uIfiotPyHzYLt9vNlg0b2R3Vh4h6Ow9dNbKNLe7cKEKv0GaMmX416ZOmo9I0RtSoNP6kT5rOmOlXd7BlCt5IxqQriR2QyrBNc4k+ZmVV31Tenf8pdru9xefau3cv+/UmDgUGMnh/MeHdlaiuk1EmTCkoKHQYdZUVzHnyt1T6dGfWVbcRVlvBjKOl3HPfHS06z/xZs/jQEEZheDRz1CouHZPWThZ7F8qEKQUFBa/HLySUcb/+Lb5HdzNhUx4lweH8TyNYt2Zzs89RXV1N/uFSdneLZtDBw1xyWWo7Wtw5UYReQUGhQ+k3/BL6j7yM/vlfkFhcx6aeSXy8fi1Hjx5t1vHZWVkURMbhUKsYaa1RwnibQBF6BQWFDmfMXffiFxzC5Wv+i49dsrxfBh98OP+c/nqXy8XGDRvZFdWbxKPHuPPmKRfI4s6FIvQKCgodjsHHl0m/fRR1TSlXrl5Fpa8/q7rH8ukni8563K5du9jpH0K10UBa4Q78Q9o+SVpXQBF6BQUFryAmKYXMK35FVOFKhu8sJT+6FysqK8jakn3GYzat30BBVC9CzXZuyex/Aa3tXChCr6Cg4DWMuPF2wmLjGbV5LlE1Dn7sm8mny76jrKzstH0rKyvJrjzGoeBQBu/MJ+0yJXb+TChCr6Cg4DVotFomPfgEwtbA1J++xKrRsTohjXnzPsJms52yb1ZWFtujeqB1uRnhqkaj03WQ1d6PIvQKCgpeRVhsPCNvnoZvyVbGZ++iMCyCTT4BLF78Ncfn/TidTtZsyWFPtzgyDxQx4YqJHWy1d6MIvYKCgtcxcPJUYgekkJr3CX3LLKzvncLaffvJzm701+/YsYNtYZHYNWqSdq4hul9iB1vs3ShCr6Cg4HUIlYqJ9z2KRqdlytqP0LkEq/sN4csl3/Dqq6/y+RdfsD2qJ/FVxxifloRQKVJ2NpT/joKCglfiHxrG2LvvR310H1eu20CJvz9ZsX2pr6+nKDicGpMvvY/sRRsd19Gmej2ajjZAQUFB4Uz0H3Ep+7I2oVq/jP59Y8iO60NBVBwWrR4h3bhUgh/Xbmb4ZaM72lSvRunRKygoeDVj77oPU1AQ4TWNpSotOgMIgRQqVvdJY1tgUAdb6P0oQq+goODVGHx9mXT/o2TH9Qdxah4bp1rDph5KJalzoQi9goKC1xM7IJU6g6HJbWdar/AzitArKCh0CoLNTSc4O9N6hZ85p9ALIQxCiE1CiK1CiO1CiP/zrO8hhNgohNgjhFgohNB51us97/d6tse3bxMUFBQuBkZuWorW6TxlndbpZOSmpR1kUeehOT16GzBGSpkKpAGXCyGGAn8DXpNSJgDHgLs9+98NHJNS9gZe8+ynoKCg0Cr67s5mwqov8K87BlLiX3eMCau+oO/uMyc9U2jknOGVsnHOcb3nrdbzksAY4BbP+jnAs8A7wFTPMsBnwD+FEEJ6Q81CBQWFTotfSCiJe/NI3Jt36vrQsA6yqPPQLB+9EEIthMgFjgLfA4VAtZTy+HNUMdDds9wdKALwbK8BQpo45wwhxBYhxJby8vLWtUJBQaHLM+qmO9Do9Kes0+j0jLqpZfVlL0aaJfRSSpeUMg2IBgYDTSV+Pt5jb6qO12m9eSnl+1LKTCllZliY8ousoKBwdvqPGs2EGQ809uCFwC80jAkzHqD/KGWy1Llo0cxYKWW1EGIVMBQIFEJoPL32aOCwZ7diIAYoFkJogACgqu1MVlBQuFjpP2q0IuznQXOibsKEEIGeZSMwDigAVgLXeXabBnzlWV7seY9n+wrFP6+goKDQcTSnRx8JzBFCqGn8YfhESrlECLED+FgI8VcgB5jl2X8WMFcIsZfGnvxN7WC3goKCgkIzaU7UTR6Q3sT6fTT663+53gpc3ybWKSgoKCi0GmVmrIKCgkIXRxF6BQUFhS6O8IZxUiFEOXDwPA8PBSra0JyORGmLd9JV2tJV2gFKW44TJ6U8Z3y6Vwh9axBCbJFSZna0HW2B0hbvpKu0pau0A5S2tBTFdaOgoKDQxVGEXkFBQaGL0xWE/v2ONqANUdrinXSVtnSVdoDSlhbR6X30CgoKCgpnpyv06BUUFBQUzoIi9ArnhRAiUAhxf0fb0dEIIS4TQizpaDuOc/J98TbbWoMQ4iEhRIEQYv4v1mcKId7sKLvaAyHEdCFEVFueUxF6hfMlELjohd4LafF98eSx8nbuByZLKW89vsKTPXeLlPKhDrSrPZgOXFxCL4TwEUJ846lZu00IcaMQ4oAQItSzPdOTOhkhxLNCiA+FEKuEEPuEEF7xARBCxAshdgoh/u1pw3whxDghxFpPzd3BQohgIcSXQog8IcQGIUSK51ivbBPwEtBLCJErhNh8cs9RCPFPIcR0z/JAIcSPQogsIcQyIURkRxn8S4QQTx3/fwohXhNCrPAsjxVCzBNCTBBCrBdCZAshPhVC+Hq2X+65n2uAazqwCU1x4r4ArwC+QojPPPbOF0IIAM936M+eNnh1biohxLtAT2CxEKJGCPG+EOI74L+d5anlDDr2Z893Z5unTUIIcR2QCcz3fLeMbWKAlNKrX8C1wAcnvQ8ADgChnveZwCrP8rPAOkBP42yzSkDrBW2IB5xAMo0/rlnAhzQWaZkKfAm8BTzj2X8MkNsJ2rTNs3wZsOSkbf+ksVei9dge5ll/I/BhR9t+kp1DgU89y6uBTR6bnwFmAj8BPp7tM4E/AwYaK6gleO7fJye3vaNfTdyXGhrrRaiA9cBIz7YDwFMdbW8L2nXA8/l/1vP9MTb12fPW1xl0LPik93OBKz3Lq4DMtry+1/fogXxgnBDib0KIUVLKmnPs/42U0ialrKCx9GG39jexWeyXUuZLKd3AdmC5bLyr+TR+OUfSeLORUq4AQoQQAZ5jvbVN56IvMAD43tPD/CONouMtZAEDhRB+gI1GIcwERgEWIBFY67F9GhAH9KPxXu7x3L95HWJ589kkpSz2fO5yafysHWdhx5jUahZLKS0dbUQLaUrHRgshNgoh8mns3CW118VbVGGqI5BS7hZCDAQmAy96Htmc/Ox2MvziENtJyy68p40n2+U+6b2bRhudpx3xcwlGb23TcU6+H/DzPRHAdinlsAtv0rmRUjqEEAeAO2l88sgDRgO9gP3A91LKm08+RgiRRhOlMb2Ys312Gi6wLW1Fp7P7DDr2Wxp77kVCiGc5XcvaDK/v0XtGn81SynnAq0AGjY9xAz27XNtBprU1PwG3QmO0BFAhpaztUIvOTh3g51k+CCQKIfSep5CxnvW7gDAhxDAAIYRWCNFuvZbz5CfgCc/f1cC9NPZ8NwAjhBC9AYQQJiFEH2An0EMI0ctz/M2nn7JDOfm+KHgJZ9AxgArP2M91J+3e5vfQ23qGTZEMvCKEcAMO4D7ACMwSQjwNbOxI49qQZ4H/CCHyADM/l2P0SqSUlZ7B5G3AtzT6qvOAPTRWHENKafcMLr3p+QHQAK/T6LryFlYDfwDWSykbhBBWYLWUstwzoPyREELv2fePnp7ZDOAbIUQFsIZG95RX8Iv7YgHKOtomBaBpHbuaRpfOAWDzSfvOBt4VQliAYW3hplJmxiooKCh0cbzedaOgoKCg0DoUoVdQUFDo4ihCr6CgoNDFUYReQUFBoYujCL2CgoJCF0cRegUFBYUujiL0CgoKCl0cRegVFBQUujj/D/8JEJGGHWsRAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# recursive multi-step forecast with linear algorithms\n",
    "from math import sqrt\n",
    "from numpy import split\n",
    "from numpy import array\n",
    "from pandas import read_csv\n",
    "from sklearn.metrics import mean_squared_error\n",
    "from matplotlib import pyplot\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.preprocessing import MinMaxScaler\n",
    "from sklearn.pipeline import Pipeline\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.linear_model import Lasso\n",
    "from sklearn.linear_model import Ridge\n",
    "from sklearn.linear_model import ElasticNet\n",
    "from sklearn.linear_model import HuberRegressor\n",
    "from sklearn.linear_model import Lars\n",
    "from sklearn.linear_model import LassoLars\n",
    "from sklearn.linear_model import PassiveAggressiveRegressor\n",
    "from sklearn.linear_model import RANSACRegressor\n",
    "from sklearn.linear_model import SGDRegressor\n",
    " \n",
    "# split a univariate dataset into train/test sets\n",
    "def split_dataset(data):\n",
    "\t# split into standard weeks\n",
    "\ttrain, test = data[1:-328], data[-328:-6]\n",
    "\t# restructure into windows of weekly data\n",
    "\ttrain = array(split(train, len(train)/7))\n",
    "\ttest = array(split(test, len(test)/7))\n",
    "\treturn train, test\n",
    " \n",
    "# evaluate one or more weekly forecasts against expected values\n",
    "def evaluate_forecasts(actual, predicted):\n",
    "\tscores = list()\n",
    "\t# calculate an RMSE score for each day\n",
    "\tfor i in range(actual.shape[1]):\n",
    "\t\t# calculate mse\n",
    "\t\tmse = mean_squared_error(actual[:, i], predicted[:, i])\n",
    "\t\t# calculate rmse\n",
    "\t\trmse = sqrt(mse)\n",
    "\t\t# store\n",
    "\t\tscores.append(rmse)\n",
    "\t# calculate overall RMSE\n",
    "\ts = 0\n",
    "\tfor row in range(actual.shape[0]):\n",
    "\t\tfor col in range(actual.shape[1]):\n",
    "\t\t\ts += (actual[row, col] - predicted[row, col])**2\n",
    "\tscore = sqrt(s / (actual.shape[0] * actual.shape[1]))\n",
    "\treturn score, scores\n",
    " \n",
    "# summarize scores\n",
    "def summarize_scores(name, score, scores):\n",
    "\ts_scores = ', '.join(['%.1f' % s for s in scores])\n",
    "\tprint('%s: [%.3f] %s' % (name, score, s_scores))\n",
    " \n",
    "# prepare a list of ml models\n",
    "def get_models(models=dict()):\n",
    "\t# linear models\n",
    "\tmodels['lr'] = LinearRegression()\n",
    "\tmodels['lasso'] = Lasso()\n",
    "\tmodels['ridge'] = Ridge()\n",
    "\tmodels['en'] = ElasticNet()\n",
    "\tmodels['huber'] = HuberRegressor()\n",
    "\tmodels['lars'] = Lars()\n",
    "\tmodels['llars'] = LassoLars()\n",
    "\tmodels['pa'] = PassiveAggressiveRegressor(max_iter=1000, tol=1e-3)\n",
    "\tmodels['ranscac'] = RANSACRegressor()\n",
    "\tmodels['sgd'] = SGDRegressor(max_iter=1000, tol=1e-3)\n",
    "\tprint('Defined %d models' % len(models))\n",
    "\treturn models\n",
    " \n",
    "# create a feature preparation pipeline for a model\n",
    "def make_pipeline(model):\n",
    "\tsteps = list()\n",
    "\t# standardization\n",
    "\tsteps.append(('standardize', StandardScaler()))\n",
    "\t# normalization\n",
    "\tsteps.append(('normalize', MinMaxScaler()))\n",
    "\t# the model\n",
    "\tsteps.append(('model', model))\n",
    "\t# create pipeline\n",
    "\tpipeline = Pipeline(steps=steps)\n",
    "\treturn pipeline\n",
    " \n",
    "# make a recursive multi-step forecast\n",
    "def forecast(model, input_x, n_input):\n",
    "\tyhat_sequence = list()\n",
    "\tinput_data = [x for x in input_x]\n",
    "\tfor j in range(7):\n",
    "\t\t# prepare the input data\n",
    "\t\tX = array(input_data[-n_input:]).reshape(1, n_input)\n",
    "\t\t# make a one-step forecast\n",
    "\t\tyhat = model.predict(X)[0]\n",
    "\t\t# add to the result\n",
    "\t\tyhat_sequence.append(yhat)\n",
    "\t\t# add the prediction to the input\n",
    "\t\tinput_data.append(yhat)\n",
    "\treturn yhat_sequence\n",
    " \n",
    "# convert windows of weekly multivariate data into a series of total power\n",
    "def to_series(data):\n",
    "\t# extract just the total power from each week\n",
    "\tseries = [week[:, 0] for week in data]\n",
    "\t# flatten into a single series\n",
    "\tseries = array(series).flatten()\n",
    "\treturn series\n",
    " \n",
    "# convert history into inputs and outputs\n",
    "def to_supervised(history, n_input):\n",
    "\t# convert history to a univariate series\n",
    "\tdata = to_series(history)\n",
    "\tX, y = list(), list()\n",
    "\tix_start = 0\n",
    "\t# step over the entire history one time step at a time\n",
    "\tfor i in range(len(data)):\n",
    "\t\t# define the end of the input sequence\n",
    "\t\tix_end = ix_start + n_input\n",
    "\t\t# ensure we have enough data for this instance\n",
    "\t\tif ix_end < len(data):\n",
    "\t\t\tX.append(data[ix_start:ix_end])\n",
    "\t\t\ty.append(data[ix_end])\n",
    "\t\t# move along one time step\n",
    "\t\tix_start += 1\n",
    "\treturn array(X), array(y)\n",
    " \n",
    "# fit a model and make a forecast\n",
    "def sklearn_predict(model, history, n_input):\n",
    "\t# prepare data\n",
    "\ttrain_x, train_y = to_supervised(history, n_input)\n",
    "\t# make pipeline\n",
    "\tpipeline = make_pipeline(model)\n",
    "\t# fit the model\n",
    "\tpipeline.fit(train_x, train_y)\n",
    "\t# predict the week, recursively\n",
    "\tyhat_sequence = forecast(pipeline, train_x[-1, :], n_input)\n",
    "\treturn yhat_sequence\n",
    " \n",
    "# evaluate a single model\n",
    "def evaluate_model(model, train, test, n_input):\n",
    "\t# history is a list of weekly data\n",
    "\thistory = [x for x in train]\n",
    "\t# walk-forward validation over each week\n",
    "\tpredictions = list()\n",
    "\tfor i in range(len(test)):\n",
    "\t\t# predict the week\n",
    "\t\tyhat_sequence = sklearn_predict(model, history, n_input)\n",
    "\t\t# store the predictions\n",
    "\t\tpredictions.append(yhat_sequence)\n",
    "\t\t# get real observation and add to history for predicting the next week\n",
    "\t\thistory.append(test[i, :])\n",
    "\tpredictions = array(predictions)\n",
    "\t# evaluate predictions days for each week\n",
    "\tscore, scores = evaluate_forecasts(test[:, :, 0], predictions)\n",
    "\treturn score, scores\n",
    " \n",
    "# load the new file\n",
    "dataset = read_csv('household_power_consumption_days.csv', header=0, infer_datetime_format=True, parse_dates=['datetime'], index_col=['datetime'])\n",
    "# split into train and test\n",
    "train, test = split_dataset(dataset.values)\n",
    "# prepare the models to evaluate\n",
    "models = get_models()\n",
    "n_input = 7\n",
    "# evaluate each model\n",
    "days = ['sun', 'mon', 'tue', 'wed', 'thr', 'fri', 'sat']\n",
    "for name, model in models.items():\n",
    "\t# evaluate and get scores\n",
    "\tscore, scores = evaluate_model(model, train, test, n_input)\n",
    "\t# summarize scores\n",
    "\tsummarize_scores(name, score, scores)\n",
    "\t# plot scores\n",
    "\tpyplot.plot(days, scores, marker='o', label=name)\n",
    "# show plot\n",
    "pyplot.legend()\n",
    "pyplot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "运行示例将评估十种线性算法并总结结果。\n",
    "\n",
    "在评估每种算法并以单行摘要报告性能时，其中包括整体RMSE以及每步RMSE。\n",
    "\n",
    "我们可以看到，大多数评估模型运行良好，整周误差均在400千瓦以下，其中随机梯度下降（SGD）回归器的表现最佳，总体均方误差（RMSE）约为383。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "还创建了针对10个分类器的每日RMSE的折线图。\n",
    "\n",
    "我们可以看到，除两种方法外，所有方法都聚集在一起，并且在整个7天的预测中表现均良好。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "通过调整某些性能更好的算法的超参数可能会获得更好的结果。此外，更新示例以测试一组非线性和集成算法可能很有趣。\n",
    "\n",
    "一个有趣的实验可能是在输入更多或更少前几天的情况下，评估一种或几种性能更好的算法的性能。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
